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1.  INTRODUCTION 

1.0  OBJECTIVES  OF  REPORT 

The  main  objectives  of  this  report  are: 

(1)  provide  the  theoretical  and  modeling  background  used  in  a 
two-dimensional  model  of  heat  and  soil -water  flow, 
coupled  by  soil-water  freezing  and  thawing  (FR0ST2X  series). 

(2)  present  the  major  modeling  assumptions  used  in  the  com¬ 
puter  program,  and  the  various  components  of  the  program 
in  order  to  aid  the  user  in  the  use  of  the  model. 

(3)  present  the  necessary  parameters  and  data  requirements 
used  in  the  computer  model. 

(4)  discuss  the  computer  model  output  product,  and  methods 
to  interpret  the  results. 

(5)  present  the  data  input  sequence  to  the  computer  model. 

1.1  REPORT  ORGANIZATION 

This  report  is  organized  into  five  chapters  and  one  appendix. 

CHAPTER  1.  Provides  an  introduction  to  the  report. 

CHAPTER  2.  Develops  the  Nodal  Domain  Integration  (NDI)  model 
of  two-dimensional  soil-water  flow  In  saturated 
and  unsaturated  soils. 

CHAPTER  3.  Develops  the  NDI  model  for  heat  flow  In  three 
dimensions.  The  third  dimension  is  Included 
In  this  NDI  development  (an  NDI  analog  for 
radial  coordinates  Is  presented  in  Appendix  A) 

In  order  to  provide  for  the  extension  of  the 
phase  change  model  to  other  coordinate  systems. 


Develops  the  two-dimensional  soil -water  phase 
change  model  used  to  couple  the  heat  and  soil- 
water  flow  during  soil-water  freezing  and  thawing. 
Presents  a  summary  of  the  two-dimensional  freezing/ 
thawing  model  and  the  modeling  input-output 
characteristics.  Provides  documentation  to  the 
computer  program  input  requirements. 

Expands  upon  the  heat  flow  modeling  approach 
developed  in  chapter  3  to  radial  coordinates. 

The  theory  developed  in  the  appendix  applies  to 
both  heat  and  soil-water  flow  in  the  selection  of 
mass- lumping  factors  for  use  in  NDI  domain  models. 

1.2  COMPUTER  PROGRAM  SOURCE  CODE 

This  report  presents  documentation  for  the  FR0ST2B  version  of  the 
FROST2X  series  of  two-dimensional  phase  change  models.  Also  provided  in 
this  report  is  the  documentation  for  the  program  PROTO0  which  provides  a 
user-friendly  data  input  capability  for  subsequent  use  by  program  FROST2B. 

Both  computer  programs  are  written  in  FORTRAN  IV  and  can  be  used  on 
most  mini-computer  class  computers. 

Computer  code  for  the  FR0ST2B  and  PROTO0  can  be  obtained  from  the  U.S. 
Army  Corps  of  Engineers,  Cold  Regions  and  Research  Laboratory  (CRREL) . 

1.3  REPORT  AUTHORIZATION 

This  report  was  prepared  under  the  direction  of  Dr.  Richard  L.  Berg 
and  Mr.  Francis  Sayles  of  CRREL,  Hanover,  New  Hampshire. 


CHAPTER  4. 


CHAPTER  5. 


APPENDIX 


2.  UNIFIED  NUMERICAL  MODEL  OF  TWO-DIMENSIONAL 
SOIL-WATER  FLOW 

2.0  INTRODUCTION 

Numeric  approximation  of  two-dimensional  non-linear  partial  differential 
equations  such  as  occur  in  the  theory  of  unsaturated  ground  water  flow  is 
generally  limited  to  numeric  solution  by  the  finite  difference  or  finite 
element  methods.  Finite  difference  approximations,  such  as  described  by 
Spalding  (1972)  for  a  specified  control  volume,  can  be  determined  for 
rectangular  and  also  for  irregular  two-dimensional  domains.  Finite  element 
methods,  such  as  variational  principle  applications  and  weighted  residuals, 
can  also  be  applied  to  irregular  two-dimensional  domains.  Both  methods 
determine  numerical  algorithms  which  are  often  compared  to  each  other  for 
numerical  "efficiency"  or  other  descriptions  of  superiority  (Hayhoe,  1978). 

Recently,  Hromadka  and  Guymon  (1982a, b)  have  developed  a  new  numerical 
approach  called  the  nodal  domain  integration  method  which  has  been  applied 
to  one  and  two-dimensional  linear  and  non-linear  problems  for  irregular 
rectangular  domains.  From  this  numerical  approach,  the  finite  difference 
and  finite  element  (Galerkin)  methods  are  "unified"  into  a  single  numerical 
statement. 

In  this  chapter  the  nodal  domain  integration  method  is  applied  to  the 
two-dimensional  triangular  finite  element.  As  special  cases,  the  finite 
element  and  finite  difference  numerical  analogs  are  determined  by  the 
appropriate  specification  of  a  single  parameter  in  the  resulting  nodal  domain 
integration  numerical  statement.  Thus,  all  three  numerical  approaches  are 
unified  by  one  numerical  statement  similar  to  the  usual  finite  element  matrix 
system. 

The  purpose  of  this  effort  is  two-fold.  The  first  objective  is  to  pre¬ 
sent  a  basic  description  of  the  nodal  domain  integration  procedure  as  applied 
to  the  class  of  partial  differential  equations  generally  encountered  in  the 
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theory  of  unsaturated  groundwater  flow.  Detailed  mathematical  derivations 
of  this  numerical  approach  are  contained  in  the  literature  (Hromadka  and 
Guymon,  1980a, b,c)  and  will  not  be  repeated  here.  The  theoretical 
foundations  of  this  numerical  method  are  based  on  the  well-known  subdomain 
technique  of  the  finite  element  weighted  residuals  approach. 

The  second  objective  is  to  present  a  unified  numerical  statement  which 
represents  the  finite  element  Galerkin  statement,  finite  difference  inte¬ 
grated  control  volume  statement,  and  the  nodal  domain  integration  statement 
(for  a  linear  shape  function  distribution  of  the  state  variable)  by  the 
specification  of  a  single  constant  parameter  in  the  resulting  triangle  element 
matrix  system. 

A  secondary  objective  is  to  briefly  discuss  the  application  of  the 
nodal  domain  integration  triangle  element  statement  to  reducing  computer 
memory  requirements  by  the  technique  of  approximating  a  higher  order  or  more 
complex  family  of  triangle  (or  global)  shape  functions  by  a  linear  shape 
function  approximation.  A  detailed  mathematical  description  of  this  linear 
shape  function  approximation  technique  is  given  for  the  one-dimensional 
case  in  a  previous  work  (Hromadka  and  Guymon,  1982a  and  will  not  be 
repeated  here. 

The  main  purpose  of  this  chapter  is  to  develop  a  generalized  unified 
triangle  element  matrix  system  which  can  be  used  in  computer  programming 
efforts.  Consequently,  a  single  computer  program  can  be  developed  which 
essentially  represents  the  finite  difference,  finite  element  (Galerkin),  and 
nodal  domain  integration  numerical  approaches.  This  work  does  not  recommend 
one  numerical  approach,  (although  for  the  problems  tested,  the  nodal  domain 
integration  scheme  produced  better  levels  of  accuracy,  i.e.  Hromadka  and 
Guymon,  1980a-c,  and  generalizes  the  numerical  expressions  so  that  the 
extension  to  the  three-dimensional  case  can  be  made  by  appropriate  inte¬ 
gration  of  the  resulting  expressions  in  the  third  dimension  (see  Chapter  3). 


2.1  GOVERNING  EQUATIONS 

Two-d imens tonal  unsaturated  Dare lan  soil-water  flow  in  a  nondefomable 
soil  matrix  fl  may  be  described  by  the  partial  differential  equation 
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where  (K  ,  K  )  are  anisotropic  hydraulic  conductivity  values  in  the  (x,y) 
x  y 

directions,  respectively;  <p  is  the  total  hydraulic  energy  head  ($  •  <|>  +  y) ; 
i|>  is  the  soil-water  pore  pressure  head;  and  6  is  the  volumetric  water 
content.  In  (1).  water  content  is  assumed  to  be  functionally  related  t>~ 
soil-water  pore  pressure  according  to  the  usual  soil  drying  curve,  with 
hysterisis  effects  neglected.  Thus, 
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where  0q  is  the  soil  porosity,  assumed  constant. 

Guymon  and  Luthln  (1974)  define  a  volumetric  water  content  to  pore  pressure 
gradient  by 
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Therefore,  (1)  may  be  rewritten  as 


(4) 
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For  ease  of  presentation,  the  soli  matrix  ft  Is  assumed  homogeneous  and 
Isotropic  with  hydraulic  conductivity  K^.  Therefore,  (4)  is  simplified 
for  discussion  purpose  as 

+57tKh!f1  -e*H<x*y)efl  <5> 

2.2  NUMERICAL  SOLUTION 

The  class  of  partial  differential  equations  including  (5)  can  be 
described  by  an  operator  relation  * 

t 

f 

A(<|>)  -  0;  (x,y)efl  (6) 

where  A(<|>)  is  the  mathematical  relation  operating  on  the  state  variable,  <J>. 

The  finite  element  method  can  be  used  to  approximate  (6)  by  the  method  of 
weighted  residuals  (Pinder  ahd  Gray,  1977). 

|  |  A($)w±dxdy  -  0  (i  «  1,2,...,M)  07) 

ft 

where  (J>  is  approximated  by  a  linear  combination  of  suitable  functions,  $, 

defined  by 
~  M 

4>  -  l  a  ♦  (8) 
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The  subdomain  version  of  weighted  residuals  uses  the  weighting  functions 
l,(x,y)e  & 
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where  ft  is  discretized  into  two-dimensional  nodal  domains  ft^ 
respect  to  finite  element  e) 
ft  -  UftJ 
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Thus,  Che  subdomain  method  is  simply  an  appropriate  integration  of  the 
governing  equations  over  each  finite  element. 

The  nodal  domain  integration  method  (Hromadlca  and  Guymon,  1980a, b) 
is  an  application  of  the  finite  element  subdomain  approach  over  appropriately 
defined  nodal  domains.  These  nodal  domains  are  partitions  of  the  finite 
element  such  that  the  resulting  nodal  domain  corresponds  to  the  control 
volume  configuration  established  by  Spalding  (1972) .  For  a  two-dimensional 
triangle  finite  element,  the  nodal  domain  partitions  are  established  by  the 
triangle  medians  as  shown  in  Fig. 2.1  allocating  one-third  of  the  triangle's 
area  to  each  respective  nodal  domain.  Detailed  mathematical  descriptions  of 
the  integration  process  for  linear  and  nonlinear  partial  differential 
operators  are  given  in  Hromadka  and  Guymon  (1980a, b,c).  These  derivations 
result  in  an  expression 
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where  Integrations  are  made  on  the  assumed  shape  function  spatial  and 
temporal  distributions  of  the  state  variable  and  nonlinear  parameters. 

This  integration  procedure  substantially  differs  from  the  integrated 
finite  difference  control  volume  approach  (Spalding,  1972)  due  to  the  finite 
difference  approach  assuming  the  state  variable  to  be  constant  interior  to 
each  control  volume  (nodal  domain)  and  due  to  the  finite  difference  approach 
assuming  flux  terms  as  constant  (spatially)  along  each  side  of  the  control 


FIG.  2.1.  LINEAR  SHAPE  FUNCTION  TRIANGULAR  FINITE  ELEMENT 
PARTITIONED  INTO  NODAL  DOMAIN  CONTRIBUTIONS 

BY  MEDIANS 


FIG.  2.2. 


VECTOR  DESCRIPTION  OF  TRIANGLE  FINITE  ELEMENT 

GEOMETRY 
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volume  rather  than  integrating  the  spatial  variation  of  the  flux  terms  as 
defined  by  the  assumed  shape  function  distribution  of  the  state  variable. 

Another  major  difference  of  the  nodal  domain  integration  method  from  the 
usual  finite  element  and  finite  difference  approaches  is  the  use  of  the 
resulting  linear  shape  function  nodal  domain  integration  numerical  statement 
to  approximate  higher  order  or  more  complex  families  of  shape  functions  with¬ 
out  an  increase  of  computer  memory  requirements  associated  to  the  use  of  such 
more  complex  families  of  shape  functions  (Hromadka  and  Guymon,  1981  ). 

In  order  to  evaluate  the  spatially  integrated  flux  terms  of  (11)  for 
each  nodal  domain  partition  of  a  triangular  finite  element,  the  triangle 
geometry  is  defined  by  a  system  of  vectors  as  shown  in  Fig. 2. 2.  For  the 
assumed  linear  shape  function  variation  of  the  state  variable,  <fr,  in  the 
finite  element  triangle  the  spatially  integrated  flux  term  contribution  to  the 
triangle  partition  associated  to  nodal  domain  ft®  is  geometrically  determined 
by  Fig. 2. 3.  From  Fig. 2.3  flux  must  contribute  to  A®  through  the  boundaries  of 
ft  |  and  can  be  calculated  by  the  flux  vector  through  state  variable  ^-values 
(at  node  1)  and  <J>*  as  shown  in  the  figure.  Thus, 

1 

♦  *  -  “  +  <(>kd2)  (12) 

The  integration  of  the  spatial  boundary  of  ft®  normal  to  the  considered  flux 
vector  is  L/2  as  shown  in  Fig. 2. 3.  Thus,  the  efflux  is  geometrically 
determined  to  be  (constant  for  linear  shape  function  approximation) 


and  th«  Integrated  efflux  contribution  for  fl®  la 


-  —  IW  +  V21  - *•  *i> 


where,  from  Fig.  2.3, 


2  -*  2  2 
L  *  r,,  •  r , ,  ■  x,.  +  y , . 
jk  jk  jk  7jk 


d,L  “  r,,  •  r -  X..X.,  +  y ... y .. 

1  ik  jk  ik  jk  'ik'jk 


d2L  ‘  rlJ  '  r)k  ■  *  (,i)xjk  +  *ij  V 


and  where  x^  *  -  x^.  In  vector  notation,  the  Integrated  efflux  contribution 

for  nodal  domain  partition  n®  of  the  triangle  finite  element  (aasumed  linear 
shape  function)  is 


<e),(  2  2, 

li)  '"jk +  V 


(xikxjk +  'ikyjk)<iij,jk +  !’iiyjk>1i 


where  Id  (18)  is  Che  area  of  Che  flnlce  elemenc  criangle,  and  Kn  ^  Che 

uniform  value  of  hydraulic  conductivity  assumed  for  a  sufficiently  small  finite 
element  triangle  (Hromadka  and  Guymon,  1980a). 

This  net  efflux  contribution  to  the  triangle  partition  of  is  identical 
to  an  Integrated  finite  difference  control  volume  approach  as  outlined  by 
Spalding.  Additionally,  the  net  efflux  term  of  (18)  is  identical  to  the 
Galerkin-determined  finite  element  triangle  component  of  the  net  efflux  for 
nodal  point  1  (Plnder  and  Gray,  1977). 


Th«  above  described  geometric  considerations  can  be  applied  for  each 
nodal  doaaln  partition  of  the  triangle  resulting  in  an  element  conduction 
matrix  identical  to  a  Galerkin-determined  element  conduction  matrix  (linear 


shape  function) 


(xjk  +  yjk*  “  (xiltXjk  +  ylkyjk)  (xijXjk  +  yijyjk^ 


Xik  +  yik 


(symmetric) 


-‘-ii-ik  +  y«ylk) 
2  .  2 
+  yU 


where  (e)  refers  to  properties  of  the  finite  element  triangle  element  (e) . 

The  notation  and  representation  of  the  element  conduction  matrix  by  (19) 
can  be  compared  to  the  finite  element  analog  given  in  Myers  (1971). 

The  nodal  domain  integration  method  evaluates  the  integration  of  the  state 
variable  in  the  triangle  partition  for  as 


$  dxdy  -  [22<0i  +  7^  +  7<J>k) 


where  (20)  represents  the  integrated  variation  of  a  linearly  distributed 
state  variable  in  (Fig.  2.4.) 

The  integrated  finite  difference  control  volume  contribution  of  the 
state  variable  integrated  or  fl®  (as  described  by  Spalding)  would  be 


$  dxdy  -  y  [ij^] 


The  Galerkln-approach  finite  element  so-called  capacitance  contribution 
is  given  by  (Myers,  1971) 


[2$.  +  +  $,,] 


In  matrix  notation,  the  nodal  dona in  integration  numerical  analog  la 
represented  by  the  element  capacitance  matrix  for  6*(e)  a  uniform  value  of  6* 
within  a  sufficiently  small  finite  element  triangle  (e) 

(23) 

where  (n  *  ^  )  corresponds  to  a  linear  shape  function  distribution  of  the 
state  variable  in  the  finite  element  triangle,  and  other  values  of  n  correspond 
to  other  shape  function  approximations  as  estimated  by  a  linear  shape  function 
approximation  (for  example,  the  alternation  theorem;  Hromadka  and  Guymon,  1981  ). 

From  (21)  and  (22),  the  element  matrix  of  (23)  also  represents  the 
Galerltin-version  of  the  finite  element  method  as  well  as  the  finite  difference 
analog  for  n  ■  (2,®),  respectively.  Thus,  the  various  numerical  approaches 
are  “unified''  by  the  nodal  domain  integration  procedure  where  the  finite  element 
and  finite  difference  approaches  are  given  by  a  specified  constant  parameter, 


,<•> 

(n) 


A(e)9*(e) 

3(n+2) 


n 

l 

l 


l  l 
n  i 
l  n 


n. 

The  resulting  nodal  domain  integration  numerical  approximation  of  (S) 
is  given  by  the  cumulative  triangle  element  contributions 


'  V 

[V 

K(e)  ■ 

+  * 
F(n) 

A. 

♦k 

where  dots  represent  time  derivatives  of  nodal  point  values  of  the  state  variable 
22 

<P,  and  n  ■  (2,  °°)  correspond  to  a  linear  shape  function  numerical  approxi¬ 

mation  by  the  Galerkin  finite  element,  nodal  domain  integration,  and  finite 
difference  methods,  respectively. 


From  the  above,  the  popular  domain  numerical  methods  of  finite  element 
and  finite  difference  are  unified  into  an  overall  numerical  approach  which  is 


a  subset  of  the  nodal  domain  integration  approach.  Consequently,  a  computer 
program  based  on  a  finite  element  approach  may  be  modified  into  the  unified 
approach  of  the  nodal  domain  integration  method. 

Nodal  domain  integration  numerical  statements  for  irregular  rectangular 
domains  (and  one-dimensional  domains)  are  contained  in  other  recent  papers 
(Hromadka  and  Guymon,  1982a, b).  Extension  to  three-dimensional  problems  is 
provided  in  Chapter  3.  Additionally  radial,  cylindrical  and  spherical  coor¬ 
dinate  systems  are  considered  in  Appendix  A. 

2.3  CONCLUSIONS 

The  nodal  domain  integration  numerical  approach  has  been  used  to  determine 
a  numerical  analog  which  incorporates  the  finite  element  (Galerkin)  and 
integrated  finite  difference  methods  as  special  cases.  The  resulting 
nunerical  statements  involve  the  same  computational  requirements  as  does  the 
finite  element  procedure.  Therefore,  the  nodal  domain  integration  procedure 
unifies  the  finite  element  and  finite  difference  approaches  by  a  single 
numerical  statement  as  a  function  of  a  single  constant  parameter.  Thus, 
computer  programs  may  be  prepared  based  on  the  nodal  domain  integration 
procedure  which  inherently  contains  both  the  finite  element  (Galerkin)  and 
finite  difference  techniques. 
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3.  UNIFIED  MODEL  OF  THREE-DIMENSIONAL  HEAT  TRAN8FER 

3.0  INTRODUCTION 


A  two-dimensional  integrated  finite  difference  method  is  presented  by 
Patankar  (1980) .  This  triangle  element  model  can  be  extended  to  a  three-dimensional 
tetrahedron  element  model,  and  both  the  two-  and  three-dimensional  finite 
element  models  can  be  compared  in  the  solution  of  linear  heat  conduction  problems. 

First,  the  integrated  finite  difference  method  is  applied  to  a  three-dimensional 
heat  conduction  process  where  the  global  domain  is  discretized  into  tetrahedra- 
shaped  finite  elements.  By  integrating  the  governing  partial  differential 
equations  on  subsets  of  the  finite  elements  (nodal  domains),  an  extension 
of  the  integrated  finite  difference  (Spalding,  1972)  analog  is  developed. 

By  using  the  subdomain  integration  version  of  the  method  of  weighted  residuals, 
another  numerical  analog  is  developed  which  is  similar  to  the  integrated 
finite  difference  approach.  Comparison  of  the  integrated  finite  difference 
and  subdomain  integration  numerical  analogs  to  one  determined  by  the  Galerkin 
finite  element  method  of  the  weighted  residuals  process  indicates  that  all 
three  analogs  determine  similar  finite  element  matrix  systems  which  when  com¬ 
bined  into  a  global  matrix  system  satisfy  both  Dirichlet  and  Neumann 
boundary  conditions. 

Hromadka  and  Guymon(1981a,b) have  examined  the  finite  difference,  sub- 
domain  integration,  and  Galerkin  finite  element  methods  for  solution  of 
partial  differential  equations  in  one-  and  two-dimensional  problems.  They 
combine  these  numerical  approaches  into  a  single  numerical  statement  which 
can  represent  any  of  the  considered  numerical  methods  by  the  specification 
of  a  single  constant  parameter,  n,  in  the  element  matrix  systems.  Examinations 
of  model  approximation  error  in  comparison  to  analytical  solutions  for  linear 
heat  and  mass  conduction  problems  indicate  that  the  element  matrix  mass 
lumping  n  must  be  variable  between  elements  and  with  respect  to  time  in 
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order  to  reduce  approximation  error.  That  is,  to  minimize  numerical  approxi¬ 
mation  error,  the  method  of  numerical  solution  generally  must  vary  through  the 
entire  range  of  domain  numerical  techniques  including  the  Galerkin  finite 
element,  subdomain  integration,  integrated  finite  difference,  (control  volume, 
or  mass  lumped  finite  element)  and  numerous  versions  of  finite  element  mass- 
weighting  schemes.  Extensive  computer  simulations  of  one-  and  two-dimensional 
problems  indicate  that  sometimes  a  single  numerical  analog  will  minimize  the 
approximation  error  for  only  a  particular  region  of  the  solution  domain  or 
only  for  a  certain  duration  of  the  simulation,  and  that  continued  use  of  the 
particular  numerical  analog  will  produce  approximation  errors  greater  than 
errors  generated  by  switching  the  method  of  numerical  solution  to  another 
technique. 

The  reduction  of  the  Galerkin  finite  element  mass  matrix  into  a  diagonal 
mass  lumped  matrix  is  well  known  (Zienkiewicz,1977) .  However,  it  is  important 
to  note  that  the  so-called  mass  lumped  diagonal  matrix  is  analogous  to  the 
integrated  finite  difference  capacitance  matrix  as  developed  by  Spalding  (1972). 

An  infinity  of  mass  weightings  of  element  nodal  contributions  can  be  determined 
directly  by  introducing  an  improved  linear  trial  function  in  the  finite  element, 
where  the  element-boundary  trial  function  continuity  requirements  are  relaxed, 
and  then  using  the  usual  subdomain  integration  version  of  the  weighted  residual 
process.  The  fact  that  certain  mass  lumping  patterns  may  improve  computational 
results  (Ramadhyani  and  Patankar,  1980)  indicates  that  a  unified  model  may  be  deve¬ 
loped  which  utilizes  a  variable  mass  matrix  and  yet  can  still  represent  the  well 
known  Galerkin,  subdomain  integration,  and  integrated  finite  difference  analogs. 

In  this  chapter,  the  general  form  of  the  three-dimensional  nodal  domain 
integration  finite  element  matrix  system  will  be  developed.  The  resulting 
element  matrix  system  will  be  shown  to  represent  an  extension  of  the  integrated 
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finite  difference  method,  subdomain  integration  version  of  the  finite  element 
method  of  weighted  residuals,  and  the  Galerkin  finite  element  method  by  the 
specification  of  a  single  constant  parameter.  The  finite  element  used  is  a 
three-dimensional  tetrahedron  with  a  linear  trial  function  used  to  approxi¬ 
mate  the  governing  transport  equation's  state  variable  within  each  element. 
Although  the  main  consideration  of  this  work  is  towards  diffusion  processes, 
an  advection  component  is  included  in  the  model  development  for  generalization 
purposes. 

The  main  objectives  in  this  chapter  are  as  follows: 

(1)  develop  a  three-dimensional  integrated  finite  difference  analog  for 
a  diffusion  (or  advection-dif fusion  process)  using  a  tetrahedron 
finite  element  discretization.  This  objective  is  an  extension  of 
the  two-dimensional  triangle  element  analog  recently  presented  by 
Baliga  and  Patankar (1980) ,  and  includes  the  development  of  a  three- 
dimensional  model  and  the  representation  of  the  model  in  finite 
element  matrix  form.  An  integrated  finite  difference  approach  is 
seen  to  result  in  a  node-centered  control  volume  analog  or,  as  it 
is  usually  referred  to  in  finite  element  models,  the  "diagonal 
mass-lumping"  finite  element  approach. 

(2)  develop  a  subdomain  integration  model  for  the  tetrahedron  finite 
element  using  a  linear  trial  function  in  each  finite  element.  The 
subdomain  integration  method  is  also  referred  to  as  a  control 
volume  approach  and  has  been  shown  to  result  in  numerical  models 
which  may  give  better  approximation  results  than  the  often  used 
Galerkin  finite  element  method.  The  subdomain  integration  analog 
will  also  be  developed  into  a  finite  element  matrix  form  in  a 
manner  such  that  the  integrated  finite  difference  and  subdomain 
integration  models  can  be  readily  compared  to  the  well  known 
Galerkin  finite  element  analog. 


(3)  develop  a  variable  "mass-lumping"  finite  element  analog. 

(4)  show  that  Dirichlet  and  Neuman  boundary  conditions  are  satisfied 
in  the  global  matrix  systems  for  each  of  the  above  determined 
analogs.  This  is  significant  due  to  the  misconception  that  the 
finite  difference  and  subdomain  integration  analogs  need  special 
formulae  to  approximate  a  Neumann  boundary  condition  (e.g.,  Bear,  1979). 

(5)  develop  a  unified  three-dimensional  domain  numerical  analog  which 
represents  all  of  the  above  numerical -models  by  the  specification 
of  a  single  parameter  in  the  resulting  unified  matrix  system. 

(6)  since  the  various  domain  approaches  determined  above  can  be  unified 
into  one  numerical  statement,  use  the  unified  model  to  examine 
model  approximation  error  by  allowing  the  unified  model  to  vary 
between  the  infinity  of  available  domain  models.  This  objective 

is  achieved  by  determining  the  optimum  capacitance  matrix  nodal  mass 
weightings  in  the  finite  element  for  two  transient  linear  heat 
conduction  problems.  Both  two-  and  three-dimensional  noual  mass 
weighting  factors  are  determined  for  comparison  purposes.  Methods 
to  determine  an  optimum  element  nodal  mass  weighting  factor,  n» 
is  the  subject  of  current  research;  however,  a  method  for  esti¬ 
mating  a  mass-lumping  factor, r|  ,  for  radial  coordinates  is  given 
in  Appendix  A,  and  in  Hromadka  and  Guymon(1981) . 

.1  GOVERNING  EQUATIONS  AND  SET  DEFINITIONS 

A  three-dimensional  advect ion-diffusion  process  in  an  inhomogeneous  aniso¬ 
tropic  nondef ormable  medium  without  sources  or  sinks  may  be  macroscopically 
described  by  the  nonlinear  partial  differential  equation 
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where  (x,y,z)  are  spatial  coordinates;  t  is  time;  (K  ■  K  ,  K  *  K ,  K  ■  K  ) 

x  XX  y  yy  z  zz 

are  principal  axis  values  of  conductivity  (e.g.  Fickian  diffusivity  or  thermal 
conductivity);  C  is  a  capacitance  coefficient;  T  is  the  state  variable  (e.g. 
temperature);  and  (U,V,W)  are  (x,y,z)-axis  advection  components  (e.g.  fluid 
velocity).  It  is  assumed  that  (1)  described  the  governing  flow  process  in 
the  nondeformable  global  domain  of  spatial  definition  £2  with  global  boundary  T. 
All  parameters  defined  in  (1)  are  assumed  variable  with  respect  to  both  space 
and  time 


?  -  C(x,y,z,t)  for  ?  e  {K,K  ,K  ,U,V,W,C)  (2) 

x  y  z 


In  vector  notation,  (1)  may  be  written  as 

f  _>  f  ” 

q  •  dr  -  C  —  dV  (3) 

J  1  3t 

r  n 

where  dT  is  the  outward  unit  normal  vector  to  surface  F,  !;  [  d  T  |  j  =  dA;  and 
3T  3T  ^  3T 

q  =  (K - UT)  i  +  (K - VT)  j  +  (K - WT)  k  (4) 

X  3x  7  3y  2  3z 

The  numerical  approximation  of  (4)  requires  a  discretization  of  the  problem 

domain,  £2.  The  subdomain  (control  volume)  and  finite  element  discretization 

processes  differ.  However  by  overlapping  the  two  main  discretization  patterns, 

the  intersection  of  subdomains  with  the  finite  elements  result  in  a  third 

discretization  composed  of  smaller  "nodal  domains"  which  are  common  to  both 

of  the  main  discretizations.  For  an  n-nodal  point  distribution  in  £2  with 

associated  subdomains  R.  and  boundaries  B, 

J  i 
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Combining  (3)  with  (5)  and  (6), 
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A  finite  element  discretization  of  ft  is  defined  by 


ft  =  Uft 


6  6 

where  each  finite  element  ft  has  a  boundary  F  where 


fte  ^  fte  =  ft®  ure 


(5) 


(6) 


(7) 

(8) 


(9) 


(10) 


(ID 


(12) 


A  set  of  nodal  domains  ft?  can  be  defined  for  each  finite  element  as  the 

J 


intersection  with  associated  subdomains 


ft®  =  R.n  ft® 


(13) 


This  set  of  nodal  domains  is  defined  for  each  finite  element  ft  by  the  index 
of  element  nodal  numbers (see  Fig.  3.1  for  the  one-diinensional  case;  the 
two-dimensional  case  is  illustrated  for  triangles  in  Fig.  3.2) 
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FIG.  3.1.  ONE  DIMENSIONAL  DISCRETIZATION  OF  GLOBAL  DOMAIN  G  SHOWING 
NODAL  POINTS  (I),  FINITE  ELEMENTS  G®,  SUBDOMAINS  Rj,  NODAL 
OOMAINS  Gj,  AND  NODES  • 


FIG.  3.2a.  FINITE  ELEMENT  ft  WITH  THREE  VERTEX 
LOCATED  NODAL  POINTS 

(3) 


(I) 

FIG.  3.2b.  FINITE  ELEMENT  PARTITIONED  INTO  NODAL  DOMAINS 


FIG.  32c.  SUBDOMAIN  Rj  AS  THE  UNION  OF  ALL  NODAL 
DOMAINS  ASSOCIATED  TO  NODAL  POINT  j. 


From  (15)  S”  is  simply  the  set  of  nodes  associated  to  element  ft*.  A  finite 
element  matrix  system  equivalent  to  the  governing  domain  equation  (1)  is 
generated  for  local  finite  element  fte  !  by 


{  f  q  •  df  -  C  —  dV  l  -  {0)  ,  j  e  S* 

1  l  e  U  3t  > 


Likewise,  a  subdomain  integration  statement  for  (1)  is  generated  for  local 
subdomain  by 


(  ^  _»  f  3T 
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l  l  l  3t 


Expanding  the  transport  integral  of  (16)  gives 
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where  (£.^,8^,8^)  are  the  direction  cosines  of  dT,  and  dA  is  the  differential 
surface  area.  Equation  (18)  can  be  rewritten  as 
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where  r6  =  boundary  of  finite  element  ft6.  The  first  integral  in  the  expansion 
of  (19)  satisfies  Neumann  boundary  conditions  on  F  or  preserves  flux  continuity 
(due  to  conduction  processes)  between  finite  elements,  ft  .  In  the  global  assem¬ 
blage  of  Uft6,  the  first  integral  in  the  expansion  of  (19)  also  satisfies  Neumann 
boundary  conditions  on  the  discretized  approximation  of  global  boundary  T  by  F  . 
From  (19),  the  element  matrix  system  of  (16)  is  given  by 


where  it  is  assumed  that  boundary  conditions  of  Neumann  or  Dirichlet  are  specified 
on  global  boundary  F. 

3.2  NUMERICAL  SOLUTIONS 

In  the  derivation  of  the  finite  element  integration  statement  of  (20)  for  ft6, 
no  specification  of  the  character  of  the  state  variable  is  assumed.  In  the  fol¬ 
lowing,  the  state  variable  T  is  assumed  to  be  adequately  approximated  by  a  linear 
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trial  function  T  in  each  finite  element  ft  .  Additionally,  each  ft  is  assumed  to 

be  a  tetrahedron  finite  element  with  four  vertex-located  nodal  points. 

Therefore 


<2n 

e  o 

where  the  L^  are  the  usual  tetrahedra  volume  local  coordinates  in  ..  ;  and 
are  nodal  point  values  of  the  trial  function  estimate  T  in  '  .  Due  to  the  linear 
definition  of  T6  in  ft6,  all  spatial  gradients  of  T6  are  constant.  Consequently, 
several  well  known  domain  numerical  solutions  of  (1)  in  ft  embodied  in  the  finite 
element  method  of  weighted  residuals  result  in  similar  numerical  approximations 
in  ft  except  for  slight  variations  in  the  element  mass  matrix. 


In  the  following,  a  Galerkin  finite  element  (linear  trial  function), 
subdomain  integration  (linear  trial  function),  and  integrated  finite  difference 
analogs  will  be  determined  and  combined  into  a  single  expression.  For  the 
finite  difference  and  subdomain  analogs,  the  integrations  of  (1)  on  a  nodal 
domain  cover  ft®  of  each  finite  element  ft*  will  be  used  to  determine  a  finite 
element  matrix  system  for  ft*.  For  all  three  numerical  methods,  the  following 
description  variable  is  defined: 
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Although  the  Galerkin  method  of  weighted  residuals  used  to  solve  (1)  in  each  ft® 
is  well  known,  its  derivation  of  a  finite  element  matrix  system  is  presented  in 
order  to  develop  some  of  the  notation  and  simplifications  used  in  the  subsequent 
determinations  of  the  subdomain  integration  and  integrated  finite  difference  analogs 

Galerkin  Method  of  Weighted  Residuals 

g 

In  local  element  ft  , 

r 

$L.dV  r  0  (2  3) 

generates  a  Galerkin  finite  element  matrix  system  for  approximation  of  (1)  on  ft6. 
Equation  (23)  is  linearized  by  assuming  all  parameters  quasi-constant  during  a 
small  timestep  At  (e.g..  Bear,  1979) 

?  =  f  +  C*  ,  C6  »  s*;  (x.y.z)  £  ft6  (24) 


For  the  linear  trial  function  T6  definition  of  the  state  variable  T  in  ft6,  the 
x-direction  terms  of  (23)  are  given  by 
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where  the  first  integral  of  (25)  satisfies  Neumann  type  boundary  conditions  on 
global  boundary  T  or  conduction  flux  continuity  between  ft  . 

For  a  linear  trial  function  T*  in  ft®  and  Neumann  boundary  conditions  on  T, 
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)  dv  (26) 
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where 


3Te  lb .  T.® 
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where  are  nodal  values  of  T  at  nodal  points  j;  V  r  volume 


of  the  tetrahedron  element;  and  b^  are  vertex  coordinate  co-factors. 


Substituting  (28)  into  (27)  gives 
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From  Figs.  3.5  and  3.4  the  shape  function  gradient  in  (29)  can  be  determined 
geometrically  as  (see  chapter  2,  equation  13  for  two-dimensions) 
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Thus, 
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3x  (hj ,x  ) 

where  (At  ,x)  is  the  projection  of  the  triangle  face  [2,3,4]  onto  the  (y,z) 
coordinate  plane. 

Simplifying  (31)  and  substituting  into  (29)  gives 
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The  (y,z)  direction  terms  are  determined  analogous  to  the  above.  The  time 
derivative  term  of  (23)  is  modeled  by 
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Solution  of  (33)  determines  the  Galerkin  finite  element  capacitance 


matrix  for  local  element  ft6 
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where  T, 


-.(j  -  (1.2, 3, 4)}  . 


Subdomain  Integration 

0  0 
A  cover  of  finite  element  Q  is  given  by  the  union  of  nodal  domains  f. 

J 

e  e 

where  j  c  S  .  For  a  tetrahedron  finite  element,  local  nodal  domain  ,.x  in  .. 

assumed  defined  by  Fig.  3.4. 


The  subdomain  integration  method  solves  (1)  in  3e  by 
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For  the  x-term  transport  components  of  (1), 
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Expanding  (36)  gives 
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The  first  integral  in  (37)  satisfies  Neumann  boundary  conditions  on  T  and  con¬ 
duction  flux  continuity  between  Qe  similar  to  the  Galerkin  formulation. 

Thus,  (36)  reduces  to 
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For  a  linear  trial  function  Te  in  Qe,  (38)  simplifies  to 
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The  surface  r*-r*^re  projection  onto  the  (y,z)  plane  is  given  by  the  integral 
with  respect  to  dA^  in  (39).  From  Fig.  3.4, 
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where(Aj  ^  is  the  projection  of  triangle  face  onto  the  (y,z)  plane. 
Additionally,  from  Fig.  3.4 
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Combining  (28),  (39),  (40)  and  (41)  gives  for  fT 
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Comparison  of  (32)  to  (42)  indicates  that  the  Galerkin  finite  element  and  sub- 

domain  integration  method  of  weighted  residuals  determine  identical  transport 
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system  matrices  for  the  assumed  linear  trial  function  T  in  .  . 

For  the  time  derivative  component  of  (36), 
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which  in  matrix  form  gives 
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Spalding  (1972)  developed  an  Integrated  finite  difference  numerical  solution 
for  partial  differential  equations  such  as  (1)  defined  on  rectangular-shaped 
control  volumes  (subdomains).  In  this  version  of  the  subdomain  integration  method, 
it  is  assumed  that 


T®  dV 


that  is,  the  nodal  value  of  T  at  node  j  is  assumed  equal  to  the  spatial  distri¬ 
bution  of  T  in  subdomain  R^. 

For  a  linear  trial  function  I  in  each  $1  ,  the  transport  terms  of  (1)  evaluated 
on  each  determines  a  conduction  matrix  identical  to  the  Galerkin  and  subdomain 
integration  approaches  previously  derived.  From  (45),  the  capacitance  matrix  from 
an  integrated  finite  difference  approach  differs  from  both  the  Galerkin  and  sub- 
domain  integration  approaches  and  is  given  in  matrix  form  by  the  so-called  lumped 


mass  diagonal  matrix 
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Nodal  Dcmain  Integration 

Hromadka  and  Guymon (1981) examined  Che  one-dimensional  form  of  (1)  and 
developed  a  modification  of  the  subdomain  integration  method  of  weighted  residuals 
which  improved  numerical  modeling  accuracy  by  the  approach  of  approximating  a 

Ag 

higher  order  trial  function  T  estimate  of  the  state  variable  T  by  using  a  linear 
trial  function  Te.  Using  this  approach,  a  significant  increase  in  modeling 
accuracy  was  achieved  while  preserving  matrix  symmetry  and  the  reduced  matrix 
sizes  associated  to  a  linear  trial  function.  By  integrating  the  governing 
equations  on  suitably  defined  nodal  domains,  a  variable  symmetric  element  capaci¬ 
tance  matrix  was  defined  which  represented  the  Galerkin  finite  element,  subdomain 
integration  and  finite  difference  methods  as  special  cases. 

0 

The  nodal  domain  integration  numerical  statement  for  solution  of  (1)  on  ft 

using  a  linear  trial  function  is  found  to  be  similar  to  the  one-dimensional  case 

0 

and  is  given  in  element  matrix  form  for  Q  by 

Ke  Te  +  Pe  (n)  T*  *  (0}  (47) 


0  0 

where  K  =  sum  of  element  f2  conduction  and  convection  matrices  given 


(for  the  x-term)  by  (32)  and  (42); 
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(48) 
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and  T  ,  T  ■  vectors  of  element  (2  nodal  values  and  time-derivative  of  nodal  values 


In  (48),  the  Galerkin  finite  element,  subdomain  integration,  and  integrated 


finite  difference  numerical  statements  for  a  linear  trial  function  in 


is  given 
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by  n  *  (2,-r^-,®),  respectively. 


Consequently,  a  single  computer  code  can  be  prepared  which  readily 
represents  each  of  the  more  popular  domain  numerical  methods  by  the  specifi¬ 
cation  of  single  parameter.  Additionally  it  should  be  noted  that  (47)  accom¬ 
modates  both  Dirichlet  and  Neumann  boundary  conditions  on  T  similar  to  the 
Galerkin  finite  element  approach 

For  the  one-dimensional  nodal  domain  integration  approach  the  corresponding 
n  parameter  was  determined  to  be  a  function  of  time  and  variable  between  finite 
elements  (Hromadka  and  Guymon,  1281b), 

n  =  n  (Qe,t)  (49) 

As  an  aid  in  selecting  the  mass-lumping  factor,  a  procedure  for  calculating 
the  factor  as  a  function  of  timestep  and  element  size  is  given  in  Appendix  A 
for  the  case  of  radial  coordinates. 

3.3  APPLICATION 

As  an  application  of  the  above  methods,  two  linear  heat  conduction 
problems  (Figs.  3.5  and  3.6)  are  modeled  to  examine  approximation  error  for 
various  values  of  mass  lumping.  For  both  problems,  a  two-  and  three-dimensional 
nodal  domain  integration  model  are  used  to  approximate  the  temperature  fields. 

Using  mean  relative  error  as  the  measurement,  the  finite  element  mass  matrix 
was  varied  by  trial  and  error  until  a  value  of  n  was  determined  such  that  the 
timestep  advancement  (Crank-Nicolson  approach)  resulted  in  a  minimum  error.  In 
all  simulations,  the  finite  element  mass-lumping  factor  (rs)  is  assumed  constant 
throughout  the  solution  domain  for  the  time  advancement  timestep.  Consequently, 
mean  relative  error  is  minimized  as  a  function  of  one  mass-lumping  variable. 

The  plots  of  n  during  the  simulation  are  given  in  Fig  3.7  where  both  two-  and 
three-dimensional  (triangle  and  tetrahedron  elements)  models  are  used  to  solve 


the  test  problems. 
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From  Fig.  3.7,  both  two-  and  three-dimensional  solutions  indicate  that  for 
the  initial  portion  of  the  numerical  solutions,  an  integrated  finite  difference 
analog  minimizes  the  error  measurement.  As  the  solution  progresses  with  time, 
however,  the  numerical  analog  approaches  a  subdomain  integration  numerical 
model.  This  variation  in  numerical  approach  was  found  similar  to  ocher  test 
problem  results  for  one-  and  two-dimensional  problems.  Apparently,  the  integrated 

finite  difference  approach  reduces  mean  relative  error  when  the  state  variable 
gradient  is  severe  within  a  finite  element,  and  a  subdomain  integration  analog 
best  serves  a  milder  variation  of  the  state  variable  in  a  finite  element. 

Use  of  different  timesteps  and  element  sizes  altered  the  shape  of  the 
n  curves  in  Fig.  3.7,  but  the  scan  of  numerical  analogs  from  integrated  finite 
difference  to  subdomain  integration  is  still  evident. 

As  suggested  by  (49),  the  n  factor  was  found  to  vary  both  between  finite 
elements  and  with  respect  to  time  for  one-dimensional  problems.  The  extension 
of  such  an  approach  to  multidimensional  problems  may  follow  based  on  the  selec¬ 
tion  of  a  rule  to  determine  n  for  each  finite  element.  One  approach  is  to 
approximate  a  higher  order  or  more  complex  trial  function  within  a  finite  element 
by  a  linear  trial  function,  and  then  use  the  improved  linear  trial  function  to 
determine  an  appropriate  n  mass-lumping  factor.  As  n  varies,  a  global  matrix 
regeneration  is  required  which  increases  computational  effort.  However  for 
highly  nonlinear  problems,  global  matrix  regeneration  may  already  be  frequently 
necessary  which  helps  offset  the  n  factor  complications.  Appendix  A  presents 
an  approach  for  selecting  the  mass-lumping  factor,  and  the  approach  tested 
for  the  case  of  a  radial  coordinate  system. 
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FIG.  3.7  OPTIMUM  MASS  WEIGHTING  7J- FACTORS  FOR  TWO  AND 
THREE-DIMENSIONAL  FINITE- ELEMENT  SOLUTION  OF  TEST 

PROBLEMS 


The  results  of  the  two  test  problems  presented  in  this  chapter  are 
typical  of  the  overall  results  obtained  by  numerous  other  computer  simu¬ 
lations  of  problems  where  analytical  solutions  exist.  Generally  speaking, 
holding  the  mass-weighting  factor  n  uniform  and  constant  in  each  finite 
element  results  in  a  numerical  method  which  provides  varying  degrees  of 
approximation  error.  That  is  for  some  problems,  a  constant  mass-weighting 
numerical  model  (such  as  Galerkin,  finite  difference,  subdomain  integration, 
or  some  other  selected  mass-weighting  model)  may  provide  a  "best"  overall 
numerical  approximation.  For  other  problems,  the  same  selected  mass-weight¬ 
ing  model  may  provide  a  "best"  numerical  approximation  for  only  a  certain 
interval  of  the  simulation  or  for  only  a  certain  region  of  the  total  problem 
domain.  This  generality  applies  to  one-,  two-  and  three-dimensional  problems 
when  using  the  well  known  rectangle,  triangle  and,  as  presented  in  this 
chapter,  tetrahedron  finite  elements. 

Although  the  test  problems  considered  here  only  indicate  the  necessary 
variation  in  uniform  finite  element  mass-weightings  to  minimize  a  measure  of 
error  (while  holding  the  mass-weighting  factor  constant  between  finite 
elements),  similar  variations  in  mass-weightings  were  found  evident  when 
holding  finite  element  mass-weighting's  constant  with  respect  to  time  and 
yet  variable  between  finite  elements.  That  is,  if  the  finite  element  mass¬ 
weighting  factor  were  allowed  to  be  variable  between  finite  elements  but  held 
constant  during  the  entire  simulation,  different  distributions  of  n  through¬ 
out  the  problem  domain  resulted  in  various  levels  of  approximation  accuracy. 

Methods  for  determining  mass-weighting  factors  for  each  finite  element 
are  presented  for  one-dimensional  problems  in  Hromadka  and  Guymon  (1981). 

Two  techniques  of  improving  a  linear  trial  function  in  each  finite  element 
are  considered.  Each  improved  trial  function  approach  results  in  a 
variable  mass-weighting  factor  for  each  finite  element. 


It  is  noted  from  the  extensive  computer  simulations  prepared  during 
this  research,  that  there  is  a  strong  correlation  between  the  minimization 
of  overall  approximation  error  and  the  use  of  a  large  n  factor  distribution 
(e.g.  a  finite  difference  analog)  in  regions  of  large  or  fast  variations 
of  the  problems?  governing  state  variable.  Additionally,  use  of  a  sub- 
domain  integration  n  distribution  in  regions  of  mild  or  slow  variations 
of  the  state  variable  also  tends  to  reduce  overall  approximation  error. 
Surprisingly,  the  various  optimized  n  curves,  such  as  shown  in  Fig.  3.7, 
did  not  include  the  often  used  Galerkin  finite  element  numerical  model. 

3.4  CONCLUSIONS 

A  nodal  domain  integration  numerical  model  is  derived  for  approximating 
a  three-dimensional  anisotropic  heat  conduction  process  in  an  inhomogeneous 
continuum.  The  nonlinear  partial  differential  equation  is  linearized  in  each 
tetrahedron-shaped  finite  element  by  assuming  all  parameters  quasi-constant 
for  small  durations  of  time.  A  linear  trial  function  is  assumed  to  adequately 
describe  the  governing  state  variable  in  each  finite  element. 

The  resulting  nodal  domain  integration  model  is  found  to  represent  the 
Galerkin  finite  element,  subdomain  integration,  and  finite  difference  methods 
as  special  cases.  Additionally,  both  Dirichlet  and  Neumann  boundary  conditions 
are  accommodated  similar  to  a  Galerkin  finite  element  numerical  model. 

Application  of  the  nodal  domain  integration  model  to  linear  heat  conduction 
problems  indicate  that  the  finite  element  mass  matrix  must  vary  with  time  in 
order  to  provide  an  optimum  numerical  solution  for  the  entire  simulation. 

The  mass  matrix  is  defined  as  a  function  of  a  single  variable,  n,  and  allows 
a  good  representation  of  the  required  mass-lumping  necessary  to  minimize 
numerical  solution  error.  Use  of  the  numerical  model  allows  a  unified  computer 
code  to  be  developed  which  offers  the  capability  to  represent  a  Galerkin,  sub- 
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domain  integration,  integrated  finite  difference,  and  an  infinity  of  different 
mass-lumped  matrix  models,  as  well  as  provide  the  capability  to  vary  between 
these  numerical  analogs  according  to  some  specified  model-selection  rule. 

From  the  complete  NDI  development  of  Chapter  3  for  three-dimensional  flow 
problems  and  the  two-dimensional  development  of  Chapter  2,  it  is  seen  that 
both  heat  and  soil-water  NDI  models  can  be  prepared  which  represent  the  Galerkin 
finite  element,  integrated  finite  difference  and  subdomain  integration 
methods  by  the  specification  of  a  single  mass-lumping  factor.  This  NDI  ap¬ 
proach  is  used  in  program  FR0ST2B  to  model  both  the  heat  and  soil-water  flow 
regimes.  Consequently,  the  program  user  can  specify  a  particular  domain  analog 
by  the  use  of  the  appropriate  mass-lumping  factor. 
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4.  ISOTHERMAL  PHASE  CHANGE  MODEL 

4.0  INTRODUCTION 

A  model  of  phase  change  in  freezing  and  thawing  soils  is  developed  for 
cold  regions  engineering  problems  which  require  two-dimensional  analysis 
of  the  thermal  regime  of  soils.  Such  problems  include  complex  boundary 
conditions  such  as  atmosphere-ground  surface  thermal  interaction  and 
snowpack-insulation.  Other  concerns  include  complex  soil  conditions  such 
as  the  presence  of  a  peaty  muskeg  or  tundra-like  soil  which  may  provide 
thermal  insulation  for  underlying  ice-rich  mineral  soil.  Although  several 
models  have  been  developed  to  predict  temperatures  in  freezing  and  thawing 
soils,  oftentimes  the  key  question  is  simply  whether  or  not  the  soil  is 
frozen  since  soil  structural  properties  are  significantly  influenced  by 
the  soil-water  state  of  phase. 

The  history  of  modeling  the  coupled  heat  and  moisture  transfer  process 
with  phase  change  in  freezing  (and  thawing)  soils  seems  to  begin  with  the 
work  of  Stefan  in  1890.  By  assuming  that  water  freezes  in  the  soil  at  a 
constant  temperature  (0°C)  the  equations  of  conduction  heat  transfer  can 
be  solved  exactly.  The  resulting  solution  is  known  as  the  moving  boundary 
problem,  and  is  inadequate  for  soils  where  moisture  is  retained  in  the 
thawed  condition  (for  freezing  temperatures),  or  when  moisture  transfer 
occurs . 

Neumann  expanded  on  Stefan’s  analysis  by  including  a  partial  differen¬ 
tial  equation  system  to  describe  the  thermal  profile  on  both  sides  of  the 
moving  boundary.  Although  inadequate  for  real  world  systems,  the  Neumann 
problem  can  be  used  to  verify  numerical  models  of  soil  freezing  wherein 
moisture  transfer  and  residual  unfrozen  moisture  contents  are  neglected. 
Such  applications  of  Neumann's  anaylsis  were  used  in  Berggren  (1943)  and 
Aldrich(1956) .  Further  exact  solutions  are  available  for  soil  freezing 
problems.  These  mathematical  developments  (Goodman,  1964;  Sikarskie  and 
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Boley,  1965)  are  valid,  however,  only  for  special  case  studies. 

For  problems  of  soil  freezing  where  an  unfrozen  moisture  content  exists 
below  the  freezing  temperatures,  a  moisture  transport  process  is  generally 
also  occuring.  Williams(1972)  concluded  that  this  moisture  transfer  effect 
is  usually  of  the  magnitude  that  ignoring  its  contribution  can  lead  to 
unacceptable  errors.  The  presence  of  moisture  below  freezing  temperatures 
negates  the  assumptions  used  in  the  moving  boundary  problems.  Thus,  an 
entirely  new  approach  to  fine-grained  soil  freezing  analysis  was  required. 
Lukianov  and  Golovko  (1957)  introduced  the  "apparent  specific  heat  capacity" 
approach  whereby  the  latent  heat  effects  of  freezing  water  are  lumped  into 
the  transient  heat  capacity  term  of  the  heat  transport  equation. 

With  the  advent  of  computer  simulation  capability,  the  past  two  decades 
have  witnessed  considerable  effort  in  freezing  soil  analysis  developments. 
These  recent  numerical  modeling  efforts  can  be  divided  into  two  groups 
depending  on  whether  soil  moisture  transfer  effects  are  included. 

Nakano  and  Brown  (1971)  modeled  the  heat  transfer  process  while 
Harlan  (1973)  approached  the  coupled  heat  and  moisture  transport  problem 
with  a  finite  difference  numerical  algorithm,  which  Is  based  on  the  assump¬ 
tion  of  an  analogy  between  moisture  transport  processes  in  unsaturated  soil 
and  a  freezing  soil.  Guymon  and  Luthin  (1974)  included  vertical  soil-water 
flow  in  a  finite  element  model  of  a  one-dimensional  vertical  soil  column. 
Sheppard  (1978)  and  Taylor  and  Luthin  (1978)  also  presented  mathematical 
models  for  the  coupled  heat  and  moisture  transfer  process.  These  models  use 
the  assumption  that  water  content  of  the  soil  is  a  function  of  temperature 
during  phase  change.  Jame  (1978)  expanded  on  Harlan's  model  to  compare 
numerically  simulated  results  to  experimental  data.  It  is  noted  that  these 
later  models  require  spatial  and  temporal  magnitudes  on  the  order  of  1.18  in, 
(3  cm)  and  0.5  hours,  respectively.  A  method  to  extend  the  above  models  to 


economically  model  large-scale,  long  duration  two-dimensional  problems  has 
not  been  advanced.  These  models  and  numerical  models  of  frost  heave  are 
reviewed  by  Guymon,  et  al.  (1980)  and  Hopke  (1980)  among  others.  It  is 
noted  that  while  one-dimensional  models  of  soil  freezing  or  thawing  are 
adequate  for  a  large  number  of  applications, at  least  two-dimensional  models 
are  required  for  many  problems  such  as  buried  pipelines,  roadway  berm  prob¬ 
lems,  and  embankments  on  permafrost. 

Geothermal  models  (i.e.,  soil-water  flow  effects  are  neglected)  often 
provide  a  thermal  analysis  capability  which  is  sufficient  for  many  problems. 
Examples  of  geothermal  models  which  are  reported  in  the  literature  include 
the  isothermal  phase  change  model  of  Bafus  and  Guymon  (1976),  the  two-  and 
three-dimensional  isothermal  phase  change  model  of  Guymon  and  Hromadka  (1977), 
and  more  recently  the  sophisticated  two-dimensional,  moving-mesh,  finite 
element  model  of  Albert  (1984). 

An  advantage  of  the  continuously  deforming  grid  system  is  that  the 
freezing  front  is  tracked  sharply,  and  requires  negligible  interpretation 
of  nodal  value  estimates  of  ice  content  to  locate  the  interface  between 
frozen  and  unfrozen  soil.  An  earlier  one-dimensional  version  of  a  finite 
element  continuously  deforming  coordinate  system  model  is  presented  in 
detail  in  O'Niell  and  Lynch  (1981).  Currently,  both  of  these  models  are 
restricted  to  homogeneous  soil  systems.  Additionally,  the  mathematical 
modeling  approach  results  in  a  type  of  convection  term  due  to  the  precise 
handling  of  the  mesh  movement.  The  mesh  moving  algorithm  requires  addi¬ 
tional  computational  effort  in  order  to  provide  a  freezing  front  coordinate- 
movement  approximation  coupled  with  an  interior  domain  nodal  point  trans¬ 
formation  and  finite  element  regeneration  capability.  Such  a  modeling 
approach  may  be  inadequate,  however,  for  soil-water  freezing/thawing 


problems  which  involve  soil-water  flow  effects  and  a  nonhomogeneous  soil 
system. 


The  main  objective  of  this  chapter  is  to  present  the  soil-water  phase 
change  modeling  approach  used  to  couple  the  heat  and  soil-water  flows 
described  in  chapters  2  and  3,  respectively.  The  ultimate  goal  of  this 
effort  is  to  develop  a  computer  model  capable  of  determining  the  thermal 
and  moisture  states  of  a  two-dimensional  soil  system  subjected  to  the 
freezing  and  thawing  processes  experienced  in  cold  climate  regions.  A 
major  factor  influencing  the  choice  of  the  modeling  approach  is  that  not 
only  should  the  computer  model  be  accurate,  it  should  also  be  economical  to 
use,  without  the  need  to  obtain  parameter  data  which  is  costly  to  evaluate. 
Additionally,  it  is  desirable  that  the  computer  code  be  prepared  such  that 
it  can  be  accomodated  on  small  computers  such  as  the  PDP  11/34  or  Data 
General  "Eclipse"  computer  systems. 


4.1  HEAT  AND  SOIL-WATER  FLOW 

The  theory  of  heat  transport  in  freezing  soils  and  their  thermal  pro¬ 
perties  have  been  the  subject  of  several  recent  publications  such  as 
Lunardini  (1981)  and  Farouki  (1981).  For  two-dimensional  heat  flow  in 
isotropic  soils,  the  governing  partial  differential  equation  is 
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where  x  and  y  =  cartesian  coordinates,  t=time,  T=soil-water-air-ice  mixture 
temperature, 9i=volume  of  ice  per  unit  total  volume  of  soil,  Ky*thermal 


conductivity  of  soil-water-air-ice  mixture,  Cm=volumetric  heat  capacity  of 
soil-water-air-ice  mixture,  L*volumetric  latent  heat  of  fusion  of  bulk  water, 
p^  -  »  density  of  ice  and  water,  respectively,  Cw=volumetric  heat  capa¬ 

city  of  water,  and  U  and  V  =the  x  and  y  velocity  flux  components.  In  (1), 
the  latent  heat  parameter,  L,  may  be  assumed  to  be  constant  for  temperatures 
less  than  -20°C  (Anderson,  et  al,  1973).  The  thermal  parameters,  K-p  and  Cm, 
are  assumed  to  be  functions  of  the  volumetric  content  of  each  material  constituent 
over  the  nodal  point  control  surface  and  control  volume  (chapter  3) .  To 
reduce  computational  effort,  the  convection  terms  in  (1)  are  approximated 
as  an  average  value  estimated  from  the  previous  timestep  solution,  and  then 
included  in  the  load  vector  term  of  the  numerical  analog.  Figure  4.1  illus¬ 
trates  the  heat  flow  conservation  model  used  for  each  nodal  point.  For 
example,  DeVries  (1966)  uses  a  volumetric  fraction  proportion 

C  -  I  C.  0  (2) 
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where  Cj=volumetric  heat  capacity  of  the  j*’*1  constituent  and  0^  =volumetric 
fraction  of  the  j**1  constituent.  The  heat  flow  equation  is  nonlinear  due 
to  the  lGj.  and  Cm  parameters  being  functions  of  ice  and  water  content  in  a 
freezing  and  thawing  soil.  Other  considerations  include  relationships 
between  temperature  and  soil-water  undergoing  phase  change  such  as  discussed 
by  Anderson,  et  al  (1973).  A  mathematical  model  of  soil-water  flow  in  un¬ 
frozen,  isotropic  soils  given  in  nondeformable  saturated  or  unsaturated 
porous  media  is  (e.g.,  Bear,  1979) 


where  $ =total  hydraulic  head,  K^=Darcy  hydraulic  conductivity. 


ALGORITHM 
STEP  NUMBER 


ALGORITHM  PROCEDURE 


Estimate  area-averaged  thermal 
parameters  of  heat  capacity  over 
the  control  volume  R . ,  and  thermal 
conductivity  on  the  dontrol 
surface  8 . . 


(2)  Estimate  nodal  values  of  tempera¬ 
ture  at  time  t+At  given  the  tem¬ 
peratures  at  time  t.  (See  Chapter 
3). 


If  nodal  temperature  values  at 
time  t  or  t+At  indicate  phase 
change  of  soil-water,  modify 
nodal  temperature  values  and 
thermal  parameters  according  to 
the  isothermal  phase  change  model 


(4)  Return  to  Step  (1)  to  model  next 
timestep  advancement. 


FIG.  4.1 .  HEAT  FLOW  MODEL 
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9U -volumetric  water  content,  and  S»moisture  sink.  Figure  4.2  illustrates  the 
soil-water  flow  conservation  model  used  for  each  nodal  point.  In  unsaturated 
soils  may  be  assumed  to  be  a  function  of  soil-water  pore  pressures  (soil- 
water  tension).  The  moisture  sink  term  for  a  freezing  soil  accounts  for 
the  phase  change  of  liquid  water  to  ice  by 
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Application  of  (3)  in  a  freezing/thawing  problem  requires  the  incorporation 
of  special  considerations  in  order  to  describe  soil-water  flow  in  a  frozen 
soil.  For  example,  the  presence  of  ice  in  soil  significantly  affects  the 
rate  of  soil-water  flow  (Nakano,  1982).  This  impact  may  be  interpreted  as 
a  reduction  in  the  Darcian  hydraulic  conductivity  (assuming  a  Darcian  soil- 
water  flow  model).  In  studies  by  Jame  (1978)  and  Taylor  and  Luthin  (1978), 
the  soil-water  flow  model  hydraulic  conductivity  had  to  be  significantly 
reduced  in  a  freezing  zone  in  order  to  adequately  reproduce  measured  thermal 
and  soil-water  data  in  freezing  horizontal  columns.  From  these  studies, 
a  general  hydraulic  conductivity  relationship  was  proposed 
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where  K=soil-water  hydraulic  conductivity,  K^=unfrozen  soil-water  hydraulic 
conductivity,  8^-volumetric  ice  content,  and  E=a  calibration  factor  evaluated 
from  one-dimensional  vertical  column  soil  freezing  laboratory  models.  In 
the  computer  model,  the  hydraulic  conductivity  in  (5)  is  assumed  to  be  zero 
for  completely  frozen  soils.  Soil-water  freezing  characteristic  curves 
indicate  that  soil-water  exists  even  at  very  cold  temperatures.  This  "resi¬ 
dual"  water  content  is  used  as  a  lower  bound  for  water  content  in  the  model. 


ALGORITHM 
STEP  NUMBER 


ALGORITHM  PROCEDURE 


(1)  Estimate  area-averaged  soil -water 
flow  parameters  of  0*  over  the 
control  volume  R, ,  hydraulic  con¬ 
ductivity  on  the'control  surface  B... 


(2)  Estimate  nodal  values  of  soil -water 
energy  head  at  time  t+At  given 
the  values  at  time  t.  (See 
Chapter  2.) 


(3)  If  phase  change  is  predicted  from 
the  HEAT  FLOW  MODEL,  modify  nodal 
values  of  water  content,  pore  pres¬ 
sure  and  energy  head  according  to 
isothermal  phase  change  model. 
Adjust  soil -water  flow  parameters 
to  accommodate  soil -water-ice 
mixture. 


(4)  Return  to  Step  (1)  to  model  next 
timestep  advancement. 
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FIG.  4.2.  SOIL-WATER  FLOW  MODEL 
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Soil-water  content  is  assumed  a  function  of  pore  pressure  head  by 


a!*|Vl 

where  0Q=porosity,  ip  =  soil -water  pore  pressure  head,  and  A  and  n  regression 
fit  coefficients. 

4.2  PHASE  CHANGE 

The  two-dimensional  phase  change  model  has  been  under  continual  develop¬ 
ment  since  1979.  The  primary  effort  has  been  under  the  direction  of  Dr. 

Richard  L.  Berg  of  the  U.S.  Army  Research  and  Engineering  Laboratory  in 
Hanover,  New  Hampshire.  The  main  objective  of  this  effort  is  to  develop 
a  computer  model  which  provides  analysis  capability  for  geotechnical  engineers 
involved  in  cold  regions  engineering  project.  The  program  is  designed  to  be 
comprehensive,  and  yet  easy  to  use  with  a  minimum  of  field  data  requirements. 
Development  of  this  model  has  been  reported  upon  as  follows  (Guymon,  et  al. 
1984):  Guymon,  et  al.  (1980)  and  Berg,  et  al.  (1980)  represented  the  con¬ 
cepts  of  a  modeling  approach  and  presented  early  verification  and  sensitivity 
results.  Subsequently,  Guymon,  et  al.  (1981a)  presented  additional  veri¬ 
fication  results,  and  Hromadka,  et  al.  (1982)  presented  a  detailed  evaluation 
of  model  sensitivity  to  the  choice  of  numerical  analog.  Guymon,  et  al.  (1981b) 
evaluate  parameter  sensitivity  and  develop  a  probabilistic  model  which  is 
cascaded  with  the  deterministic  one-dimensional  model.  Finally,  Guymon,  et  al. 
(1984)  present  applications  of  the  two-dimensional  model  applied  to  pipeline 
studies  and  roadway  embankments. 

A  summary  of  the  key  assumptions  used  in  the  two-dimensional  model  are 


as  follows: 


1.  Soil-water  flow  occurs  in  unfrozen  zones  by  liquid  water  film 
driven  by  hydraulic  gradients  such  as  described  by  Darcy's  law. 

2.  Soil-water  flow  in  frozen  zones  is  negligible. 

3.  Heat  flow  is  mathematically  described  by  (1). 

4.  Soil-water  flow  is  mathematically  described  by  (3). 

5.  Soil-water  phase  change  may  be  approximated  as  an  isothermal 
process . 

6.  Unfrozen  zones  are  nondeformable,  and  in  freezing  zones  or  frozen 
zones,  deformation  is  due  to  ice  segregation  or  lens  thawing  only. 

7.  Soil  water  pore  pressures  in  freezing  zones  are  governed  by  a 
residual  water  content  determined  from  soil  freezing  tests. 

8.  Hysteresis  effects  are  negligible. 

9.  Salt  exclusion  processes  are  negligible. 

10.  Constant  parameters  (e.g.,  porosity)  remain  constant  with  respect 
to  time;  i.e.,  freeze-thaw  cycles  do  not  modify  parameters. 

11.  Freezing  and  thawing  processes  in  a  two-dimensional  medium  occur 
in  such  a  way  that  there  are  no  internal  shear  or  stress  forces 
developed  between  different  zones. 

The  phase  change  model  which  couples  the  heat  and  soil-water  flow  models 
is  based  upon  the  simple  control  volume  approach  where  freezing/thawing  occurs 
isothermally  (Guymon  and  Hromadka,  1977).  The  algorithm  is  based  on  a  control 
volume  generated  by  the  nodal  integration  method.  A  volume  of  freezing  soil 
is  retained  at  0°C  until  the  latent  heat  of  fusion  of  all  available  soil-water 
for  freezing  has  been  evolved. 

Figure  4.3  illustrates  the  phase  change  model  logic,  and  indicates  how 
the  entire  nodal  point  control  volume  is  lumped  during  phase  change. 
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FIG.  4.3  ISOTHERMAL  SOIL-WATER  FREEZING  MODEL 


I 
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4.3  MODEL  PARAMETER  RELATIONSHIPS 


Because  the  soil-water  flow  model  is  based  on  energy  head,  solution  of 
(3)  requires  that  the  relationship  between  soil-water  pressure  head  and 
soil-water  content  be  defined  where  the  time-derivative  term  is  replaced  by 

30.  _  36  9^  (7) 

3t  It 

39 

The  partial  derivative,  is  determined  from  Gardner's  (1958)  relationship 

of  (6) .  Hydraulic  conductivity  is  defined  as  a  function  of  by 


Aj*r+i 

where  KQ=saturated  hydraulic  conductivity,  and  Ag  and  m  are  calibrated 
parameters.  Hydraulic  conductivity  in  soil-freezing  zones  is  estimated  by 
the  exponential  relationship  of  (5) . 

Soil-water  pore  pressures  at  freezing  fronts  are  defined  by 

ip  (at  freezing  front)  =  i p  (8rES) 

where  0res  is  a  residual  soil-water  content. 

Additional  data  required  for  use  with  the  model  are  boundary  conditions 
for  the  heat  and  soil-water  flow  models,  and  initial  conditions  of  temperature, 
water  content  and  ice  content.  Boundary  conditions  are  modeled  as  constants 
or  sinusoidal  functions  specified  by  the  program  user  as  to  amplitude  and 
wave  period.  Modeling  data  input  requirements  are  discussed  in  detail  in 
Chapter  5, 
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6.  PROGRAM  PROTON 

5.0  INTRODUCTION 

A  FORTRAN  computer  program  is  available  which  accomodates  two-dimen¬ 
sional  heat  and  soil-water  flow  models  (Chapters  2  and  3)  as  coupled  by 
an  isothermal  phase  change  model  (Chapter  4).  The  program  can  be  used  to 
analyze  two-dimensional  freezing/thawing  problems  which  have  sufficient 
known  information  to  supply  the  necessary  modeling  parameters,  boundary 
conditions,  and  initial  conditions. 

Because  of  the  sophistication  of  the  two-dimensional  phase  change 
model  and  the  data  requirements  needed  to  properly  represent  inhomogeneity 
of  the  system,  boundary  conditions,  and  other  complexities,  a  special  data 
input  program  is  developed  in  order  to  aid  the  model  user.  This  general 
purpose  data  preparation  program,  PROTO0,  develops  the  data  input  file  to 
be  used  directly  by  the  two-dimensional  phase  change  program. 

In  this  chapter,  the  PROTO0  program  will  be  reviewed  in  detail.  The 
actual  CRT  screen  pages  will  be  displayed  which  show  the  data  entry  prompts 
in  their  respective  order  of  appearance.  Also  shown  on  the  CRT  screen 
pages  are  the  computer  program  variable  names  associated  to  each  prompt  in 
order  to  aid  in  understanding  the  FORTRAN  code. 

5.1  MODELING  APPROACH 

Heat  and  Soil-Water  Flow 

Chapters  2  and  3  provide  the  details  of  the  mathematical  models  used 
to  approximate  the  thermal  and  soil-water  effects  in  two-dimensional  problems. 
The  analyst  initially  discretizes  the  problem  geometry  (two-dimensional 
domain)  into  a  collection  of  triangle  finite  elements  (Fig.  5.1).  The 
computer  model  further  subdivides  these  triangles  into  nodal  domains  such 
as  shown  in  Fig.  5.2.  The  collection  of  nodal  domains  forms  a  control  volume 
for  each  nodal  point  (Fig.  5.3).  The  computer  model  then  balances  the 


heat  and  soil-water  flow  over  each  nodal  control  volume  using  straight-line 
interpolation  of  temperature  and  soil-water  energy  head  to  compute  the 
corresponding  rates  of  flow.  This  straight-line  interpolation  function  is 
shown  graphically  in  Fig.  5.4. 

Phase  Change 

Chapter  4  describes  the  isothermal  phase  change  algorithm  used  to 
approximate  the  freezing  and  thawing  of  soil-water.  This  algorithm  is 
based  upon  the  "lumped-mass"  control  volume  shown  in  Fig.  5.5  (which 
conforms  geometrically  to  the  control  volume  used  to  balance  heat  and  soil- 
water  flow).  Figure  5.6  illustrates  the  budget  used  for  each  nodal  control 
volume  which  accounts  for  the  residual  water  content  (unavailable  for 
freezing),  the  remaining  water  content  available  for  freezing,  and  the  ice 
content.  During  phase  change,  the  nodal  temperature  is  defined  to  equal 
the  freezing  point  depression  (usually  0°C) . 

5.2  MODELING  PARAMETERS 

The  parameters  needed  for  the  computer  model  fall  into  3  categories: 

(i)  heat  flow  parameters 

(ii)  soil -water  flow  parameters 

(iii)  phase  change  parameters 
Heat  Flow  Parameters 

Thermal  conductivities  and  heat  capacities  are  required  to  model  heat 
flow.  These  parameters  can  be  usually  developed  from  published  formulas, 
or  obtained  from  charts  and  tables. 

Soil-Water  Flow  Parameters 

The  soil-water  flow  model  requires  information  regarding  the  coefficient 
and  exponent  used  in  the  Gardner's  function  (see  chapter  4 )  relationships 
for  hydraulic  conductivity  and  water  content  as  related  to  pore  pressure. 
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FIG.  5.5.  LUMPED  MASS  CONTROL  VOLUME 
(USED  FOR  ICE  CONTENT  BALANCE) 


FIG.  5.6.  WATER  CONTENT  BUDGET  FOR  NODE  j 


LEGEND 
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-  HEAT  AND  SOIL -WATER 
FLOW  ACROSS  BOUNDARY 
(CONDUCTION  PARAMETERS 
USED) 


-  CONTROL  VOLUME  AREA 
(CAPACITANCE  PARAMETERS 
USED) 


FIG.  5.7.  FLOW  BALANCE  MODEL  AND  PARAMETER  USAGE 


Phase  Change  Parameters 

The  volumetric  latent  heat  of  fusion  and  freezing  point  depression  are 
defined  by  the  program  user. 

Parameter  Groups 

Rather  than  enter  the  several  parameters  for  each  finite  element, 
parameter  groupings  are  used  to  combine  similar  properties.  Two  types  of 
groupings  are  used;  namely,  element  parameter  groups  and  nodal  parameter 
groups . 

-Element  parameter  groups  include  the  conduction  parameter 
data.  This  enables  a  better  description  of  conduction 
values  used  to  compute  flow  rates  across  control  volume 
boundaries. 

-Nodal  parameter  groups  include  the  capacitance  parameter 
data.  This  enables  the  program  user  to  best  define  those 
parameters  which  are  averaged  over  the  control  volume 
area  (e.g.,  heat  capacity). 

Figure  5.7  demonstrates  the  flow  balance  models  used,  and  where 
the  conduction  and  capacitance  parameters  are  assumed  to  apply. 


5.3  MODELING  RESULTS 
Nodal  Values 

The  computed  results  provide  nodal  point  values  of  temperature,  volu¬ 
metric  water  content,  and  volumetric  ice  content  produced  at  time  intervals 
specified  by  the  program  user.  A  special  feature  afforded  by  the  program  is 
the  ability  to  also  print  the  previous  timestep  computed  results  (along  with 
the  current  modeling  results)  in  order  to  compare  the  change  in  nodal  values 
of  several  variables  during  the  recent  timestep  advancement. 


Freezing  Front  Interpretation 

Because  the  ice  content  values  are  specified  at  nodal  points,  and  due  to 
the  mass-lumping  budget  used  for  the  phase  change  algorithm  (see  Fig,  5.6), 
interpretation  of  the  nodal  ice  content  values  are  required  in  order  to  locate 
the  freezing  front  (i.e.,  the  line  separating  the  frozen  soil  from  unfrozen 
soil).  Similarly,  the  temperature  values  require  interpretation  in  order  to 
locate  the  0°C  isotherm  (freezing  front)  at  the  boundary  between  the  frozen 
and  unfrozen  soil  regions.  The  interpretation  effort  required  is  directly 
related  to  the  size  of  the  finite  elements  used.  Large  triangular  elements 
necessarily  result  in  large  control  volume  dimensions  associated  to  each 
nodal  point  (see  Figs.  5.2  and  5.3). 

The  usual  interpretation  procedure  is  to  simply  assume  that  the  volume 
of  ice  estimated  to  exist  in  a  nodal  control  volume  exists  as  a  single  piece, 
and  is  located  to  the  side  of  the  control  volume  which  has  frozen. 

This  interpretation  can  be  illustrated  in  terms  of  a  one-dimensional 
problem  involving  rectangular-shaped  finite  elements.  Figure  5.8a  shows 
a  nodal  point  control  volume  which  is  initiating  freezing  of  available  soil- 
water.  The  most  recent  timestep  only  evolved  enough  heat  to  freeze  10- 
percent  of  the  soil-water  available  for  phase  change.  Due  to  the  lumped- 
mass  model,  the  entire  nodal  control  volume  is  associated  with  the  nodal 
point  value;  hence,  the  nodal  values  of  temperature  and  ice  content  indicate 
the  freezing  point  depression  and  10-percent  frozen  soil-water,  respectively. 
Thus  this  information  must  be  interpreted  to  indicate  that  the  10-percent 
frozen  available  soil-water  is  in  one-piece  and  located  as  shown  in  Fig.  5.8a. 
Figure  5.8b  illustrates  the  same  control  volume  with  60-percent  of  the 
available  soil-water  frozen.  Figure  5.8c  illustrates  a  two-dimensional 
control  volume  with  60-percent  of  the  available  soil-water  frozen. 
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FIG.  5.8c.  60- PERCENT  FROZEN  CONTROL  VOLUME  (TRIANGLE  ELEMENTS) 


5.4  THE  TWO-DIMENSIONAL,  PHASE  CHANGE  PROGRAM  SYSTEM 

The  two-dimensional,  coupled  heat  and  soil-water  flow,  with  an  isothermal 
phase  change  approximation  computer  model  is  available  as  the  FR0ST2X  series 
of  programs  where  currently  four  versions  are  available  (Table  5.1). 

Similarly,  the  PROTO0  program  has  been  extended  to  serve  special  purpose 
problems  where  the  problem  geometry  is  of  a  typical  character.  For  example, 
PR0T01  enables  for  a  quick  data  file  preparation  for  roadway  embankment 
problems  where  the  interior  nodal  points  and  finite  elements  are  developed 
by  the  program  based  on  the  entry  of  a  few  critical  geometric  coordinates  of 
the  problem  boundary  and  locations  of  regional  homogeneity  (i.e.,  identical 
parameters  for  heat  and  soil-water  flow). 

It  is  noted  that  by  considering  the  data  entry  requirements  used  in 
PROTO0,  the  engineer  can  prepare  special  purpose  data  file  preparation  codes 
which  are  compatible  with  the  FR0ST2X  series,  significantly  reducing  the 
data  entry  requirements  associated  with  the  PROTO0  general  purpose  code. 

5.5  PROTO0  DATA  REQUIREMENTS 

The  data  entry  requirements  associated  with  PROTO0  fall  into  four  broad 
categories.  These  data  groupings  are  illustrated  in  Fig.  5.9. 


TABLE  5.1  PROGRAM  FROST2X  DESCRIPTIONS 


PROGRAM 

FROST2A 

FROST2B 


FROST2C 


DESCRIPTION 

Base  Development  Version  of  FROST2X  Series. 

Two-Dimensional  Heat  and  Soil-Water  Flow 
Model  With  Isothermal  Phase  Change  Model. 
Accomodates  Heterogeneous  But  Isotropic 
Soil  Systems.  Includes  an  Apparent  Heat 
Capacity  During  Phase  Change  Compatible 
With  PROTO0  For  Data  File  Preparation. 

Extends  FR0ST2B  To  Include  Anisotropic  Soil- 
Water  Flow. 


FR0ST2D  Extends  FR0ST2C  To  Include  Vertical  Frost 

Heave  and  Overburden  Effects.  Compatible 
With  PR0T01  For  Data  File  Preparation  of 
Roadway  Embankment  Problems 


PROGRAM 
FROST 2 B 


FROST 2D 


TABLE  5.2  SUBROUTINE  TABULATION 

Main,  Indata,  Bcsine,  Trans,  Phase, 
Output,  Presol,  Comb,  Finsol. 

Main,  Indata,  Bcsine,  Trans,  Phase, 
Output,  Presol,  Comb,  Finsol,  Over. 


Note:  The  PROTO0  program  presented  in  this  chapter  is  compatible  with 
FR0ST2B. 


•  MODEL  CONTROL  DATA 

nodal  domain  integration  mass  lumping  factor 
timestep 

time  between  global  matrix  regeneration 
time  of  simulation 
time  between  output  of  results 
model  selection: 

heat  and  soil-water  flow 
heat  flow 
soil-water  flow 

include  isothermal  phase  change 
thermal  parameters  of  water  and  ice 
number  of  nodes 

number  of  temperature  boundary  condition  nodes 
number  of  pore  pressure  boundary  condition  nodes 
number  of  triangle  finite  elements 

•  FINITE  ELEMENT  PARAMETER  GROUPS 

thermal  conductivity  of  soil 
heat  capacity  of  soil 

saturated  hydraulic  conductivity  of  soil 

exponent  in  Gardner's  hydraulic  conductivity 

coefficient  in  Gardner's  hydraulic  conductivity 

hydraulic  conductivity  ice  content  correction  factor  (exponent) 

•  NODAL  POINT  PARAMETER  GROUPS 

soil  porosity 

exponent  in  Gardner's  water  content 
coefficient  in  Gardner's  water  content 
frozen  soil  residual  water  content 
heat  capacity  of  control  volume 

•  FINITE  ELEMENT  NODE  NUMBERS  AND  (x,y)  COORDINATE  DATA 

NODAL  POINT  INITIAL  CONDITIONS 
temperature 

soil-water  pore  pressure  head 
ice  content 

•  TEMPERATURE  BOUNDARY  CONDITIONS 

node  number 
maximum  temperature 
minimum  temperature 
sine  period 
phase  shift 

•  SOIL-WATER  PORE  PRESSURE  BOUNDARY  CONDITIONS 

node  number 

maximum  pore  pressure  head 
minimum  pore  pressure  head 
sine  period 
phase  shift 


FIG.  S.9.  PROGRAM  FR0ST2B  DATA  REQUIREMENTS 


S.6  PROTO0  DATA  ENTRY 


Program  PROTO0  prompts  the  model  user  for  all  data  entries.  In  the 
following,  the  data  entry  prompts  are  shown  in  their  order  of  appearance. 
Included  with  the  prompts  are  the  associated  PROTO0  variable  names.  It 
is  noted  that  several  of  the  prompts  include  suggested  parameter  values 
(for  typical  soil-water  phase  change  problems),  and  the  range  of  values 
allowed  for  use  with  program  FR0ST2B. 


*************************************************************************** 

TWO-DIMENSIONAL 

HEAT  AND  SOIL  WATER  FLOW  MODEL 
WITH  AND  WITHOUT  PHASE  CHANGE 


DEVELOPED  AT 

UNIVERSITY  OF  CALIFORNIA*  IRVINE. 


PRINCIPAL  INVESTIGATOR!  GARY  L.  GUYMON 
PROGRAM  DEVELOPMENT!  TED  V.  HROMADKA 


VERSION  DATE!  APRIL*  1982 

till  ENTER  A  C1D  TO  CONT I NUE ••«•»••• 1 


**************************ttt*******: ***********************t*tttt*********t 

♦PROGRAM  PROCESS  SELECTION* 

**********************************tMtt************************$$t********** 


ENTER  PROGRAM  PROCESS  NUMBER! 

1  =  CREATE  A  NEW  DATA  BANK 

2  =  CONTINUE  CREATING  A  DATA  BANK 

3  =  EDIT  AN  EXISTING  DATA  BANK 

4  =  EXIT  PROGRAM  CPROTOOA 3 


KOPT 


iLi’l 


introduction: <protooa> 

THIS  GENERAL  PURPOSE  DATA  FILE  PREPARATION  PROGRAM 

requires: 

1 • NO  MORE  THAN  90  NODAL  POINTS. 

2.  NO  MORE  THAN  ISO  TRIANGULAR  FINITE  ELEMENTS. 
(THREE  VERTEX-LOCATED  NODES  TO  EACH  ELEMENT) 

3.  ALL  COORDINATES  ARE  IN  FIRST  QUADRANT. 

4.  ALL  NODAL  POINTS  ARE  NUMBERED  SEQUENTIALLY 
FROM  NUMBER  1. 

5.  ALL  ELEMENTS  ARE  NUMBERED  SEQUENTIALLY  FROM 
NUMBER  1. 

6.  MATRIX  BANDWIDTH  MUST  NOT  EXCEED  IB 
(BANDWIDTH  =  MAXIMUM  ARITHMETIC  DIFFERENCE  BETWEEN 

ANY  FINITE  ELEMENT  NODAL  NUMBERS.  ♦  1) 


:::enter  a  ci3  to  continue: : : : : 


ENTER  NODAL  DOMAIN  INTEGRAT ION( NDI )  MASS-LUMPING  FACTOR: 

note:mass-lumping  factors  for  numerical  methods  are: 

2  «  GALERKIN 

3  *  SUBDOMAIN  INTEGRATION 

1000  *  INTEGRATED  FINITE  DIFFERENCE 
(ALLOWABLE  VALUES  ARE  BETWEEN  Cl  3  AND  C 100000  3) 


XNETA 


ENTER  TIMESTEP  MAGNITUDE  (HOURS): 

(CRANK-NICOLSON  TIME  ADVANCEMENT  METHOD  IS  USED  FOR 
HEAT  TRANSPORT  MODEL.  FULL Y- IMPL 1C  I T  TIME  ADVANCEMENT 
METHOD  IS  USED  FOR  SOIL-MOISTURE  TRANSPORT  MODEL) 


DELT 


::::::::::::::::::::::: 


DEFINITION  OF  'UPDATE*  i 

IN  THIS  COMPUTER  MODEL  ALL  THERMAL  AND  SOIL-MOISTURE 
PARAMETERS  ARE  HELD  CONSTANT  FOR  A  USER-SPEC  IF IFD 
DURATION  OF  TIME.  THIS  DURATION  OF  TIME  IS  CALLED  THE 
THE  'UPDATE  PERIOD*  AND  IS  EXPRESSFD  IN  UNITS  OF  HOURS 


ENTER  UPDATE  PERIOD  (HOURS) : 

(SMALLEST  ALLOWABLE  VALUE  IS  t 

: : : : : : : :::::: : : : : : : : : : : : : : : : : : : :  i 


0.1000003 ) 


UPDAT 


THE  TOTAL  COMPUTER  MODEL  SIMULATION  TIME  IS 
EXPRESSED  IN  UNITS  OF  DAYS 


ENTER  LENGTH  OF  SIMULATION  (DAYS).* 

(SMALLEST  ALLOWABLE  VALUE  IS  C  0. 0041673 ) 


EN8IM 


COMPUTER  SOLUTIONS  ARE  PRINTED  ACCORDING  TO  A  USER 
SPECIFIED  SIMULATION  INTERVAL . THIS  OUTPUT  TIME 
INTERVAL  IS  EXPRESSED  IN  UNITS  OF  DAYS 


ENTER  TIME  INTERVAL  BETWEEN  COMPUTER  OUTPUTS  (DAYS): 
(SMALLEST  ALLOWABLE  VALUE  IS  C  0.0041671) 
(LARGEST  ALLOWABLE  VALUE  IS  C  2.00001) 


OUT 


NDI  MASS  LUMPING  FACTOR  a  1000.000 
TIME  STEP  MAGNITUDE  ■  0.100000  HOURS 

UPDATE  PERIOD  =  0.100000  HOURS 

LENGTH  OF  SIMULATION  =*  2.000000  HAYS 

TIME  INTERVAL  BETWEEN  COMPUTER  OUTPUTS  = 


(INPUT  EXAMPLE) 


1.000000  DAYS 


: : select  Data  option  number.'  : : : : : : : : : : : : : : : : : 

1  =  ACCEPT  DATA  AND  CONUNUE 

2  =  DATA  IS  UNACCEPTABLE*  REINPUT  LAST  SFQUENCE 

3  =  ACCEPT  DATA  AND  TERMINATE  PROCESS 


*************************************************************************** 

PROGRAM  MODEL  SELECTIONS 

*************************************************************************** 


ENTER  PROGRAM  MODEL  OPTION  NUMBER.* 

0  *  HEAT  AND  MOISTURE  TRANSPORT 

1  *  HEAT  TRANSPORT  ONLY 

2  =  MOISTURE  TRANSPORT  ONLY 


NPATHI 


IN  THIS  MODEL •  SOIL  WATER  PHASE  CHANGE  IS  ASSUMED 
TO  OCCUR  ISOTHERMALLY.  SOIL-WATER  AND  ICE  CONTENTS 
ARE  ACCOUNTED  FOR  BY  A  SIMPLE  CONTROL  VOLUME 
BUDGET  KEEPING.  A  CONTROL  VOLUME  IS  DEFINED  FOR 
EACH  NODAL  POINT  BY  THE  SUMMATION  OF  EACH 
TRIANGLE  ELEMENT  NODAL  DOMAIN  ASSOCIATED  TO 
EACH  NODE. 


ENTER  FROGRAM  SOIL-MOISTURE  FREEZING/THAWING  MODEL 
OPTION  NUMBER.* 

0  *  EXCLUDE  ISOTHERMAL  FREEZING/THAWING  MODEL. 
1  =  INCLUDE  ISOTHERMAL  FREEZING/THAWING  MODEL. 


IPHASE 
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USER  SPECIFIED  MODEL  SELECTIONS 
HEAT  AND  MOISTURE  TRANSPORT 
INCLUDE  FREEZING/THAWING  MODEL 


(INPUT  EXAMPLE) 


•  • *  SELECT  OAT  A  OPT I ON  NUMBER  *•*•••••{•»•»»««•»•»«{••«•*•••••••!»••«•• 

1  *  ACCEPT  DATA  AND  CONTINUE 

2  *  DATA  IS  UNACCEPTABLE i  REINPUT  LAST  SEQUENCE 

3  *  ACCEPT  DATA  AND  TERMINATE  PROCESS 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

♦SPECIFIED  THERMAL  INFORMATION# 


sxxxxxxs 


ENTER  FREEZING  POINT  DEPRESSION  OF  WATER  (C  DEGREES):  TFPD 

note: 

ENTER  COI  FOR  METRIC  SYSTEM 
(ALLOWABLE  VALUES  ARE  BETWEEN  t-1003  AND  CIOOI) 


VOLUMETRIC  HEAT  CAPACITY  OP  WATER J < cal/c»**3 ) 
NOTE J ENTER  Cl. 003  FOR  METRIC  SYSTEM 


c 


: 


ENTER  VOLUMETRIC  HEAT  CAPACITY  OF  ICE !< cal /c»**3 ) 
NOTEJENTER  CO. 463  FOR  METRIC  SYSTEM 


•  • 


•  • 


•  •  •  • 


ENTER  THERMAL  CONDUCTIVITY  OF  WATER J < cal /hr . ca . c >  JK 

NOTE  I  ENTER  C4.83  FOR  METRIC  SYSTEM 


ENTER  THERMAL  CONDUC T IV I T i  OF  ICE J < c a  1 /hr . c* . c >  T| 

NOTEJENTER  C19.03  FOR  METRIC  SYSTEM 


ENTER  VOLUMETRIC  LATENT  HEAT  OF  F US  I  ON :  ( c a i /r»* *  3 >  ) 

NOTEJENTER  C80.J  FOR  METRIC  SYSTEM 


f  ^SPECIFIED  THERMAL  INFORMATION* 

THE  THERMAL  FREEZING  POINT  DEPRESSION  OF  WATER  = 
HEAT  CAPACITY  OF  WATER  =  1.0000 

HEAT  CAPACITY  OF  ICE  =  0.4400 

THERMAL  CONDUCTIVITY  OF  WATER  =  4.8000 

THERMAL  CONDUCTIVITY  OF  ICE  =  19.0000 
LATENT  HEAT  OF  FUSION  =  80.000000 


(INPUT  EXAMPLE) 


0.000000 


.'SELECT  DATA  OPTION  NUMBER .' .' .'  .*  J  .*  .* .' .'  i 

1  =  ACCEPT  DATA  AND  CONTINUE 

2  =  DATA  IS  unacceptable;  RETNPUT  LAST  SEQUENCE 

3  =  ACCEPT  DATA  AND  TERMINATE  PROCESS 


NODAL  DOMAIN  INTEGRATION 
MODEL  PROBLEM-DEFINITION  INFORMATION 


*****************************************«***M*************************** 
♦ELEMENT  PARAMETER  GROUPING  INFORMATION* 
*************************************************************************** 


IN  THIS  COMPUTER  MODEL  t FINITE  ELEMENT  PARAMETER 
INFORMATION  IS  GROUPED  INTO  'PARAMETER  GROUPINGS'  WHEREBY 
THE  APPROPRIATE  GROUPING  NUMBER  IS  USED  TO  SPECIFY 
PARAMETER  INFORMATION  FOR  EACH  FINITE  ELEMENT . THEREFORE » 
THE  USER  SPECIFIES  A  GROUPING  NUMBER  RATHER  THAN  ENTERING 
OFTEN-REPETATIUE  PARAMETER  DATA  FOR  EACH  ELEMENT.  UP  TO 
C 103  ELEMENT-PARAMETER  GROUPINGS  MAY  BE  DEFINED. 


ENTER  THE  NUMBER  OF  ELEMENT-PARAMETER  GROUPINGS? 
(ALLOWABLE  VALUES  ARE  BETWEEN  C13  AND  C103> 


NEPQ 


ENTER  ELEMENT-PARAMETER  GROUPING  *C  1  3 1 NFORMAT ION  S 

ENTER 

THERMAL  CONDUCTIVITY  OF  SOIL  <  ca  1 /hr  .  c*. .  c  >  ! 

PARELE(J.I) 

ENTER 

VOLUMETRIC  HEAT  CAPACITY  OF  SOIL  ( ca 1 /cm**3 ) ? 

P  ARELE(  J,2) 

ENTER 

THE  SATURATED  HYDRAULIC  CONDUCTIVITY  OF  SOIL  <  cm/hr ) ? 

P  ARELE(  J,3) 
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fUlin  ItHHIUMUtf 


SIMILAR  TO  THE  'ELEMENT  PARAMETER  GROUPINGS' » FURTHER 
PARAMETER  INFORMATION  IS  DEFINED  FOR  NODAL  POINTS. 

THIS  'NODAL  POINT  PARAMETER  GROUPING'  INFORMATION  IS 
IDENTIFIED  BY  A  NODAL -PARAMETER  GROUPING  NUMBER. 

THESE  SETS  OF  PARAMETERS  ARE  USED  TO  MORE  PRECISELY 
DEFINE  ASSUMED  CONSTANT  PARAMETERS  ON  A  CONTROL 
VOLUME  BASIS  FOR  EACH  NODAL  POINT  RATHER  THAN 
SIMPLY  AVERAGING  FINITE  ELEMENT  PARAMETERS 
FOR  EACH  NODAL  CONTROL  VOLUME. 

UP  TO  CIO]  NODAL-PARAMETER  GROUPINGS  MAY  BE  DEFINED. 


ENTER  THE  NUMBER  OF  NODAL-PARAMETER  GROUPINGS i  NNPG 

(ALLOWABLE  VALUES  ARE  BETWEEN  Cl]  AND  CIO]) 


ENTER  NODAL-PARAMETER  GROUPING  #C  1  ] INFORMATION 5 


•  •  ♦ 


•  •  •  • 


ENTER  THE  SATURATED  VOLUMETRIC  MOISTURE  CONTENT  OF  SOIL!  PARNOD(J.l) 
(ALLOWABLE  VALUES  ARE  BETWEEN  CO]  AND  C1J  cm**3/cin** 3 ) 


ENTER  MULTIPLIER  OF  PORE  PRESSURE  HEAD  IN  GARDNERS 
VOLUMETRIC  MOISTURE  CONTENT  FUNCTION! 


P  ARNOD(  J,  2) 


P  ARNOD(  J,3) 


SOIL  WATER  IS  AVAILABLE  FOR  PHASE  CHANGE.  THIS 
UNAVAILABLE  SOIL  WATER  IS  CALLED  A  VOLUMETRIC 
UNFROZEN  WATER  CONTENT  FACTOR. 


ENTER  THE  ASSUMED  VOLUMETRIC  UNFROZEN  WATER  CONTENT  PARNOD(J,4) 

FACTOR  OF  THE  SOIL? 

(ALLOWABLE  VALUES  ARE  BETWEEN  COT  AND  C1T) 


ENTER  THE  VOLUMETRIC  HEAT  CAPACITY  OF  THE  SOIL  <cal/c»**3 )  .'  PARNOD(J.G) 


NODAL-PARAMETER  GROUPING  *  1  (INPUT  EXAMPLE) 

SATURATED  VOLUMETRIC  MOISTURE  CONTENT  OF  SOIL  =  0.4000 

MULTIPLER  OF  PORE  PRESSURE  IN  GARDNERS  MOISTURE  CONTENT  FUNCTION  =  0.0040 

EXPONENT  OF  PORE  PRESSURE  IN  GARDNERS  MOISTURE  CONTENT  FUNCTION  =  1 . 2000 

VOLUMETRIC  UNFROZEN  WATER  CONTENT  FACTOR  =  0.1500 

VOLUMETRIC  HEAT  CAPACITY  OF  THE  SOIL  -  0.5000 


: :  .'select  data  option  number:  J  J  ;  i  j  : :  s : : : : : : : : : : : : : : 

1  *  ACCEPT  DATA  AND  CONTINUE 

2  *  DATA  IS  UNACCEPTABLE »  REINPUT  LAST  SEQUENCE 

3  =  ACCEPT  DATA  AND  TERMINATE  PROCESS 


*************************************************************************** 

•PROGRAM  PROCESS  SELECTION* 

(PROGRAM  CPROTOOB3 ) 

a************************************************************************** 


ENTER  PROGRAM  PROCESS  NUMBER.*  KOPT 

1  ■  CREATE  A  NEW  DATA  BANK  (CONTINUATION  OF  PROTOOA) 

2  *  CONTINUE  CREATING  A  DATA  BANK  (CONTINUATION  OF  PROTOOB) 

3  -  EDIT  AN  EXISTING  DATA  BANK 

4  =  EXIT  PROGRAM  CPROTOOB3 


A****.********************************************************************* 

•FINITE  ELEMENT  INFORMATION* 

*************************************************************************** 


THE  PROBLEM  DOMAIN  IS  ASSUMED  DISCRETIZED  INTO 
TRIANGLE-SHAPED  FINITE  ELEMENTS.  SINCE  THIS 
DATA  PREPARATION  PROGRAM  IS  FOR  AN  ARBITRARY 
DOMAIN.  THE  USER  MUST  ENTER  ALL  GEOMETRIC 
INFORMATION  INCLUDING  THE  PROBLEM  MESH. 

EACH  FINITE  ELEMENT  HAS  3  NODAL  POINTS  .THESE 
VERTEX  NODAL  POINTS  MUST  BE  LISTED  IN  COUNTER 
CLOCKWISE  ORDER. 


.*.*:  TENTER  A  C13  TO  CONTINUE 


ENTER  ELEMENT  NUMBER  C  13  INFORMATION 


ENTER  THE  NODE  NUMBER  OF  THE  FIRST  VERTEX  1 


IDTELE(J.I) 


ENTER  THE  NODE  NUMBER  OF  THE  SECOND  VERTEX:  IDTELE(J,2) 


ENTER  THE  NODE  NUMBER  OF  THE  THIRD  VERTEX: 


IDTELE(J,3) 


ENTER  ELEMENT-PARAMETER  GROUP  NUMBER:  IDTELE(J,4) 


ELEMENT  NUMBER  II 


(INPUT  EXAMPLE) 


NODE  NUMBER  OF  THE  FIRST  VERTEX  = 
NODE  NUMBER  OF  THE  SECOND  VERTEX  = 
NODE  NUMBER  OF  THE  THIRD  VERTEX  *  ; 

ELEMENT-PARAMETER  GROUP  NUMBER  -  1 


: : : select  data  option  number: : ::::::::::::::::::::: 

1  *  ACCEPT  DATA  AND  CONTINUE 

2  =  DATA  IS  UNACCEPTABLE*  REINPUT  LAST  SEQUENCE 

3  =  ACCEPT  DATA  AND  TERMINATE  PROCESS 


*************************************************************************** 

ENTER  INITIAL  TEMPERATURE  .INITIAL  PORE  PRESSURE  HEAD. 

AND  INITIAL  VOLUMETRIC  ICE  CONTENT  OF  EACH  NODE  AS 

requested: 

*************************************************************************** 

ENTER  INITIAL  TEMPERATURE  AT  NODE  .*  1  TOLD(J) 

«  t  t  •  t  «•  t  t  t  ••  I  «•  t  •«  t  t  »«•  I  t  •••  I  t  •  4  |  »  4  •  t  •••»•««  t  «  t  ««««•••  |  t  .  « 

ENTER  INITIAL  PORE  PRESSURE  HEAD  AT  NODE  :  1  POLD(J) 

ENTER  INITIAL  VOLUMETRIC  ICE  CONTENT  AT  NODE  S  1  XICEOL(J) 


NODE  NUMBER  1 : 


(INPUT  EXAMPLE) 


INITIAL  TEMPERATURE  OF  THE  NODE  =  l . OOOO 
INITIAL  PRESSURE  HEAD  OF  THE  NODE  =  -10.0000 

INITIAL  VOLUMETRIC  ICE  CONTENT  OF  THE  NODE  =  0.0000 


:::select  data  option  number::::::::::::::::::::::: 

1  =  ACCEPT  DATA  AND  CONTINUE 

2  *  DATA  IS  UNACCEPTABLE?  REINPUT  LAST  SEQUENCE 

3  =  ACCEPT  DATA  AND  TERMINATE  PROCESS 


♦SPECIFIED  BOUNDARY  CONDITION  INFORMATION* 


NOTE J ALL  SPECIFIED  BOUNDARY  CONDITION  < TEMPERATURE  AND 
PORE-PRESSURE  HEAD)  ARE  EXPRESSED  IN  TERMS  OF  A 
SINUSOIDAL  VARIATION. 


ENTER  NODE  NUMBER  WITH  SPECIFIED  THERMAL  BOUNDARY  CONDITION  NBCTtJ) 
<  AND  SPECIFY  THE  TEMPERATURE  FUNCTION): 


ENTER  MAXIMUM  TEMPERATURE  ON  THE  SINE  CURVE  DESCRIBING 
THE  BOUNDARY  CONDITION  AT  NODE  !  1 


BCT(J.I) 


ENTER  MINIMUM  TEMPERATURE  (IN  THE  SINE  CURVE  DESCRIBING  BCT(J  2) 

THE  BOUNDARY  CONDITION  AT  NODE  J  1  ' 


ENTER  THE  PERIOD(HOURS)  OF  THE  SINE  CURVE  DESCRIBING  BCT(J  3) 

THE  TEMPERATURE  BOUNDARY  CONDITION  AT  NODE  :  1  ’ 


ENTER  THE  PHASE  SHIFT(HOURS)  OF  THE  SINE  CURVE. 
DESCRIBING  THE  TEMPERATURE  BOUNDARY  CONDITION  AT  NODE  : 


BCT(  J,4) 


THREE  POSSIBLE  VALUES  CAN  BE  USED  TO  DESCRIBED  THE 
SHIFTING  SINE  CURVE  ! 

1.  ENTERING  A  ‘O'  VALUE:  THIS  IS  A  STANDARD  SINE  CURVE. 

2.  ENTERING  A  ‘POSITIVE*  VALUE!  THIS  SINE  CURVE  IS  BEING  SHIFTED  BACKWARD . 
BY  THAT  ENTERED  VALUE »  FROM  THE  STANDARD  SINE  CURVE. 

3.  ENTERING  A  ‘NEGATIVE*  VALUE!  THIS  SINE  CURVE  IS  BEING  SHIFTED  FORWARD. 
BY  THAT  ENTERED  VALUE,  f ROM  THE  STANDARD  SINE  CURVE. 


NODE  NUMBER  WITH  SPECIFIED  TEMPERATURE 


1 


(INPUT  EXAMPLE) 


MAXIMUM  TEMPERATURE  ON  THE  SINE  CURVE  *  -5.0000 
MINIMUM  TEMPERATURE  ON  THE  SINE  CURVE  ■  -5.0000 
THE  PERIOD  OF  THE  SINE  CURVE  -  48.0000 
THE  PHASE  SHIFT  OF  THE  SINE  CURVE  *  0.0000 


:::select  data  option  number: ::::::::::::::::::::::::::: : ::::: : :::::: : 

1  -  ACCEPT  DATA  AND  CONTINUE 

2  =  DATA  IS  UNACCEPTABLE?  REINPUT  LAST  SEQUENCE 

3  -  ACCEPT  DATA  AND  TERMINATE  PROCESS 


ENTER  NODE  NUMBER  WITH  SPECIFIED  SOIL-MOISTURE  NBCP(J) 

BOUNDARY  CONDITION  (AND  DEFINE  THE  PORE-PRESSURE  HEAD 
FUNCTION  AT  THAT  NODAL  POINT)! 

*  * . .  t  *  »  *  ♦  •  •«#•  i  I  *****♦••*••♦•««*♦*••♦**#•*  I  ♦  t  »♦  I  I  l  , 


ENTER  MAXIMUM  PORE  PRESSURE  HEAD  ON  THE  SINE  CURVE  BCP(  J,  1 ) 

DESCRIBING  THE  BOUNDARY  CONDITION  AT  NODE  :  41 


ENTER  MINIMUM  PORE  PRESSURE  HEAD  ON  THE  SINE  CURVE  BCP(J,2) 

DESCRIBING  THE  BOUNDARY  CONDITION  AT  NODE  !  41 


ENTER  THE  PERIOD ( HOURS )  OP  THE  SINE  CURVE  DESCRIBING 
THE  PORE  PRESSURE  HEAD  BOUNDARY  CONDITION  AT  NODE  S  41 


BCPCJ.3) 


ENTER  THE  PHASE  SHIFT(HOURS)  OP  THE  SINE  CURVE 

DESCRIBING  THE  PORE  PRESSURE  HEAD  BOUNDARY  CONDITION  AT  NODE  :  41 


THREE  POSSIBLE  VALUES  CAN  BE  USED  TO  DESCRIBED  THE  BCP(J,4) 

SHIFTING  SINE  CURVE  : 

1.  ENTERING  A  *0*  VALUE.*  THIS  IS  A  STANDARD  SINE  CURVE. 

2.  ENTERING  A  'POSITIVE'  VALUES  THIS  SINE  CURVE  IS  BEING  SHIFTED  BACKWARD, 
BY  THAT  ENTERED  VALUE,  FROM  THE  STANDARD  SINE  CURVE. 

3.  ENTERING  A  'NEGATIVE*  VALUES  THIS  SINE  CURVE  IS  BEING  SHIFTED  FORWARD, 
BY  THAT  ENTERED  VALUE,  FROM  THE  STANDARD  SINE  CURVE. 


NODE  NUMBER  WITH  SPECIFIED  PORE  PRESSURE  HEAD  =  41  (INPUT  EXAMPLE) 

MAXIMUM  PORE  PRESSURE  HEAD  ON  THE  SINE  CURVE  =-10.0000 
MINIMUM  PORE  PRESSURE  HEAD  ON  VHF  SINE  CURVE  =-10.0000 
THE  PERIOD  OF  THE  SINE  CURVE  =  48.0000 

THE  PHASE  SHIFT  OF  THE  SINE  CURVE  =  0.0000 


S  S  SELECT  DATA  OPTION  NUMBERS  SSSISSSSSSSSSSSSSSSSSS! 

1  =  ACCEPT  DATA  AND  CONTINUE 

2  =  DATA  IS  UNACCEPTABLE;  REINPUT  LAST  SEQUENCE 

3  =  ACCEPT  DATA  AND  TERMINATE  PROCESS 


i 
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5.7  APPLICATIONS 


| 

[  Three  example  problems  are  presented  which  illustrate  the  data  file 

development  by  use  of  PROTO0,  and  the  performance  of  the  program  FR0ST2B. 
Heat  Flow  Application 

A  one-dimensional  domain  of  unit  length  is  discretized  by  8  triangle 
finite  elements  as  shown  in  Fig.  5.10a.  At  time  t=0,  the  temperature  is 
given  by  T(x)=l.  Boundary  conditions  are  given  at  x=0  and  x=l  by 
T(x=0)=T(x=l)=0.  Using  a  normalized  timestep  of  t=0.01,  the  computed 
results  from  FR0ST2B  and  the  exact  solution  (Myers,  1971)  are  shown  in 
Fig.  5.10b.  The  data  file  prepared  by  PROTO0  is  shown  in  Fig.  5.11. 
Soil-Water  Flow  Application 

A  vertical  homogeneous  soil  column  is  discretized  by  triangle  finite 
elements  as  shown  in  Fig.  5.12a.  A  water  table  forms  the  base  of  a  steady 
state  45-degree  pore-pressure  head  profile  through  the  vertical  column. 

The  column  is  insulated  on  both  sides.  The  top  of  the  column  is  suddenly 
flooded  at  a  uniform  depth  of  2  cm.  of  water.  The  FROST2B  modeling  results 
are  shown  in  Fig.  5.12b.  The  data  file  developed  from  PROTO0  is  shown  in 
Fig.  5.13. 

Phase  Change  Model  Application 

The  vertical  column  of  Fig.  5.12a  is  now  considered  with  respect  to 
soil-  water  freezing.  Initially,  the  column  is  at  a  uniform  temperature  of 
+0.1°C.  The  top  of  the  column  is  suddenly  set  at  a  constant  temperature  of 
-5°C,  The  FR0ST2B  modeling  results  are  shown  in  Fig.  5.14.  The  data  file 
developed  by  PROTO0  is  shown  in  Fig.  5.15. 
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FIG.  5.10a  HEAT  FLOW  EXAMPLE  PROBLEM 


FIG.  5.10b  HEAT  FLOW  EXAMPLE  PROBLEM  SOLUTION 
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FIG.  5.11.  PROTO 0  DATA  FILE  FOR  HEAT  FLOW  PROBLEM  (lot  2) 
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FIG.  5.12a  SOIL- WATER  FLOW  EXAMPLE  PROBLEM 
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FIG.  5.12b  SOIL- WATER  FLOW  EXAMPLE  PROBLEM  SOLUTION 
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FIG.  5.13.  PROTO0  DATA  FILE  FOR  SOIL -WATER  FLOW  PROBLEM 

(lot  3) 
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FIG.  5.15.  PROTO 0  DATA  FILE  FOR  PHASE  CHANGE  PROBLEM  (lot  3) 


14 

13 

16 

1 

14 

11 

12 

1 

12 

15 

14 

1 

13 

14 

17 

1 

17 

16 

13 

1 

17 

14 

15 

1 

15 

18 

17 

1 

16 

17 

20 

1 

20 

19 

16 

1 

20 

17 

18 

1 

18 

21 

20 

1 

19 

20 

23 

1 

23 

22 

19 

1 

23 

20 

21 

1 

21 

24 

23 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0.10000 

0.00000 

0.00000 

0. 10000 

0.00000 

0.00000 

0. 10000 

0.00000 

0.00000 

0.10000 

-5.00000 

0.00000 

0.10000 

-5.00000 

0.00000 

0.10000 

-5.00000 

0.00000 

0.10000 

-10.00000 

0.00000 

0.10000 

-10,00000 

0.00000 

0.10000 

-1 0.00000 

0.00000 

0.10000 

-15.00000 

0.00000 

0.10000 

-15.00000 

0.00000 

0. 10000 

-15.00000 

0.00000 

0.10000 

-20,00000 

0.00000 

0.10000 

-20.00000 

0.00000 

0.10000 

-20,00000 

0.00000 

0.10000 

-25.00000 

0.00000 

0.10000 

-25.00000 

0.00000 

0.10000 

-25.00000 

0.00000 

0.10000 

-30.00000 

0.00000 

0.10000 

-30.00000 

0.00000 

0.10000 

-30.00000 

0.00000 

0.10000 

-35,00000 

0.00000 

0.10000 

-35.00000 

0.00000 

0.10000 

-35.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0,00000 

0.00000 

1 

0.1000 

0.1000 

48.0000 

0.00000 

2 

0.1000 

0.1000 

48,0000 

0.00000 

3 

0. 1000 

0.1000 

48.0000 

0.00000 

22 

-5.0000 

-5.0000 

48.0000 

0.00000 

23 

-5.0000 

-5.0000 

48,0000 

0.00000 

24 

-5.0000 

-5.0000 

48.0000 

0.00000 

0 

0.0000 

0.0000 

0.0000 

0.00000 

0 

0.0000 

0.0000 

0.0000 

0.00000 

1 

0.0000 

0.0000 

48.0000 

0.00000 

2 

0.0000 

0.0000 

48.0000 

0.00000 

3 

0.0000 

0.0000 

48.0000 

0.00000 

0 

0.0000 

0.0000 

0.0000 

0.00000 

FIG. 5.15.  (2  of  3) 


101 


0*0000  0.0000  0.0000 
111 
111 
0  0  1 

28  24  6 

0  0  0 

0  0  0 

0  0  24 


0.00000 


FIG.  5.15.  (3  of  3) 


APPENDIX  A:  NDI  Model  For  Radial  Coordinates 


A.O  INTRODUCTION 

In  chapter  3  Is  presented  the  development  of  a  nodal  domain  integration 
(NOI)  model  of  three-dimensional  heat  conduction  based  on  a  tetrahedron  finite 
elenent.  From  this  numerical  model,  the  finite  difference,  subdomain  Integra¬ 
tion  and  Galerkln  finite  element  methods,  and  an  infinity  of  finite  element 
mass-lumped  matrix  models  are  unified  Into  a  single  numerical  statement. 

As  an  extension,  the  NDI  technique  is  now  applied  to  the 
radial -coordinate,  finite  element.  It  is  shown  that  the  Galerkln  finite 
element,  subdomain  Integration,  and  an  integrated  finite-difference  numerical 
models  are  obtained  by  the  appropriate  specification  of  a  single  parameter  in 
the  resulting  NDI  statement.  .  Thus,  all  three  numerical  approaches  are  unified 
into  one  numerical  statement  similar  in  form  to  a  Galerkln  finite  element 
matrix  system.  The  extension  of  the  NDI  technique  to  developing  unified  cylin¬ 
drical  and  spherical  coordinate  models  follows  from  the  derived  radial  coordinate 

model  and  the  NDI  three-dimensional  tetrahedron  finite  element  model  of 
chapter  3, 


The  purpose  of  this  section  is  three-fold.  The  first  objective  is  to  present 

a  basic  description  of  the  NDI  technique  as  applied  to  the  class  of  partial 
differential  equations  generally  encountered  in  the  theory  of  diffusion  and 
heat  conduction.  Sufficient  set  definitions  and  integral  manipulations  are 
provided  in  order  that  the  extensions  of  the  results  to  cylindrical  and  spherical 
coordinate  systems  are  direct.  The  theoretical  foundations  of  this  numerical 
technique  are  based  on  the  subdomain  version  of  the  finite  element  weighted 
residuals  approach,  and  incorporate  the  mass-lumping  techniques  used  in  some 
finite  element  approaches.  The  development  closely  follows  the  presentation  of 
Chapters  2  and  3. 
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The  second  objective  is  to  develop  the  NOI  numerical  statement  which 
represents  the  finite  element  Galerkln  statement,  subdomain  integration 
numerical  statement,  and  integrated  finite-difference  control  volume  statement, 
by  the  specification  of  a  single  parameter  in  the  resulting  radial  element 
matrix  system. 

A  third  objective  is  to  use  the  unified  NDI  formulation  to  gain  Insight 
as  to  the  performance  of  the  several  well-known  domain  models  in  the  approxi¬ 
mation  of  radially  symmetric  heat  conduction  processes.  Since  the  NDI  model 
represents  each  of  the  most  popular  numerical  models  as  point  values  of  the 
NOI  approximation  statement,  the  same  computer  code  can  be  used  for  each 
numerical  analog,  as  well  as  an  arbitrary  finite  element  mass-lumping  scheme. 

A.l  NDI  MODEL  DEVELOPMENT 

The  partial  differential  equation  describing  radially  symnetric  heat 
conduction  in  an  isotropic  homogeneous  medium  is  given  by 

3  ["  9$  9<t 

—  Irk  —  =  rs  —  (i) 

3R  |_  3R_l  9t 

where  K  is  the  thermal  conductivity;  S  is  the  heat  capacity;  R  is  the  radial 
coordinate;  <fi  is  temperature;  and  t  is  time. 

The  finite  element  technique  apnroximately  solves  the  governing  equation 
on  a  finite  element  discretization  of  the  domain  (Zienkiewicz,  1977).  The 
integrated  finite-difference  method  uses  a  control  volume  discretization 
(Spalding,  1972). 

The  nodal  domain  integration  approach  partitions  a  finite  element  into 
smaller  "nodal  domains"  which  are  geometrically  defined  as  the  intersection 
of  the  finite  elements  and  control  volumes.  The  utility  of  this  further 
partitioning  of  the  finite  element  is  that  an  integrated  finite-difference  or 
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a  subdomain  Integration  analog  can  be  conveniently  written  In  terms  of  a  matrix 
system  similar  to  the  Galerkln  finite  element  matrix  system.  Additionally, 
flux  type  boundary  conditions  can  be  accomodated  on  the  problem  domain 
boundary  without  the  need  of  special  equations  or  finite-difference  approximations. 
The  following  set  definitions  of  subdomains  ( control - vol umes ) ,  finite  elements, 
and  nodal  domains  will  be  used  to  develop  the  NDI  finite  element  matrix  system. 
Consider  the  partial  differential  operator  relationship 

A($)  *  f;  (x,y)  e  ft,  n»nUr  (2) 

defined  on  global  domain  ft  with  boundary  condition  types  of  Dlrlchlet  or 
Neumann  specified  on  global  boundary  r.  A  n-nodal  point  distribution  can 
be  defined  In  ft  with  arbitrary  density  such  that  an  approximation  9  for  $ 
is  defined  in  ft  by 

n 

i  *  l  Mx.y)#.;  (x,y)  eft  (3) 

j»l  J  J 

where  Nj(x.y)  are  linearly  independent  global  shape  functions  and  <j>j  are 
assumed  values  of  the  state  variable,  <j>,  at  nodal  point  j.  In  (3)  it  is 
assumed  that  except  for  a  set  of  Lebesque  measure  zero 

lim  $  *  lim  $  =  <t>,  (x,y)  eft  (4) 

n-»  max  1 1  ( x  ■  ,y  ■ )  ,  ( \  »yk )  1 1 

A  closed  connected  spatial  subset  R.  is  defined  for  each  nodal  point  j  such 

*j 

that 

n  =  U  **  (5) 

j-1  J 
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The  sets  Rj  are  generally  called  control -volumes  or  subdomains,  and  usually 
are  accoapanled  with  additional  requirements  that 

(xj,  yj)  e Rj;  (xj,  yj)  i  Rk,  jfk  (6) 


and 


(7) 


where  (xy  y^)  are  the  spatial  coordinates  of  node  j  and  Bj  Is  the  boundary 
of  Rj.  It  Is  assumed  that  every  subdomain  is  disjoint  except  along  shared 
boundaries,  l.e. 

RjnRk-BjfiBk  (s) 

The  subdomain  method  of  the  finite  element  weighted  residuals  approach 
approximates  (2)  by  solving  the  n  equations 


where 


|  { A( <t>)  -  f)  Wj  dA  »  0,  j  *  l,2,***,n  (9) 

n 


1.  (x.y)  e  Rj 
.  0,  (x.y)  {  Rj 


(10) 


A  second  cover  of  Q  is  defined  by  the  finite  element  method  with 

n  *  Une  (li) 


where  ne  is  the  closure  of  finite  element  domain  Qe  and  its  boundary  re. 


Let  St  b«  tha  set  of  subscripts  defined  by 

Se  5  {JlflTlRj  A  {*»  (12) 

that  Is,  Se  Is  simply  the  set  of  nodes  associated  to  ft*. 

Then  a  set  of  nodal  domains  ft*  Is  defined  for  each  finite  element  domain  ft* 
by  (Fig.  A.l) 

n*  »  fle  0  Rj .  J  e  Se  (13) 

The  subdomain  method  of  weighted  residuals  as  expressed  by  (9)  can  be 
rewritten  In  terms  of  the  subdomain  cover  of  ft  by 


’ 

(AU)  -  f)  Wj  dA 

ft 


(AU)  -  f)  dA 


With  respect  to  the  finite  element  discretization  of  ft. 


(14) 


(AU)  -  f)  dA  *  (AU)  -  f)  dA  (15) 

(  •  R .  Oft* 

J  J 


where  for  each  finite  element  ft*  a  matrix  system  is  given  by  generatinq 
for  each  nodal  point  j  eSg 


(AU)  -  f)  dA 


ft*OR 

J 


(A($)  -  f)  dA, 


J 

ft 


e 

j 


(16) 


From  the  above  subset  definitions  and  set  covers  of  ft,  application  of  the 
usual  subdomain  method  to  the  governing  partial  differential  operation  of 
(2)  is  accomplished  by  an  Integration  of  the  governing  equations  over  the 
nodal  domains  interior  of  each  finite  element,  resulting  in  a  finite  element 
matrix  system  similar  to  that  determined  by  the  Galerkin  finite  element  method. 
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A. 2  NDI  INTEGRATION 

In  this  section,  the  governing  heat  flow  equation  Is  Integrated  over  the 
several  nodal  doaalns  associated  to  a  finite  element  ne.  This  approach  Is 
simply  the  subdomain  Integration  weighted  residual  method  as  applied  to  a  sub- 
domain  or  control  volume,  except  that  the  approximation  error  Is  averaged 
over  the  nodal  domains  Interior  of  the  subdomain.  These  nodal  domain  contri¬ 
butions  can  then  be  reassembled  Into  matrix  form  for  each  element  fte.  Using 
the  previous  set  notation,  the  operator  relationship  for  the  radial ly-synmetrlc 
heat  conduction  model  of  (1)  Is 


A4)-f 


a 

3R 


(17) 


whereby  substituting  (17)  Into  (16)  gives  an  element  matrix  system  for  fte 


f  a 

3R 


l  p 

Expanding  (18), 


3$ 
RK  — 
3R 


a*  ) 

RS  — 

at 


dA 


*  (0)  ,jeSe 


(18) 


’ 

3$ 

30 

RK  — 

ds 

.  +  . 

RK  — 

ds 

i 

an 

j 

3n  , 

r.t 

!Dre 

r 

e 

Le 

.per 

Ire 

r 

i 

J  J 


30 

RS  —  dA 
3t 


.  J  t 


n' 


(19) 


where  the  first  term  of  (19)  cancels  due  to  flux  contributions  from  con¬ 
tiguous  finite  elements  or  satisfies  Neumann  boundary  conditions  on  global 
boundary  r ;  and  where  (n,s)  are  normal  and  tangential  vector  components  on 
Bj ,  r*.  and  r*. 


\  AA  VjVjS  .■  *•  .N  \  ,  *.  *.  *, 
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In  order  to  evaluate  the  Integration  expressions  of  (19),  definitions 
of  the  finite  element  end  subdomain  covers  of  global  domain  ft  are  necessary. 
The  finite  element  cover  ft*  of  ft  Is  assumed  defined  by  the 


ft1  =  ((x,y)|0  ■  r,  <  R  <  rf) 
ft2  =  {(x,y)|r2  <  R  <  r2} 

ft""1  s  {(x.y)|rn-1  <  R  <  rn  -  L) 
where  r.  Is  the  radial  coordinate  of  node  1;  and 

ft  =  { ( x ,y ) 1 0  <  R  <  L). 


(20) 


The  subdomain  cover  Rj  of  ft  Is  assumed  defined  by 

R,  =  {(x,y)|0  «  2ri  <  2R  <  {r;  ♦  r2 ) } 

Rz  =  ((x,y)|(rl  ♦  r2)  <  2R  <  (p^  ♦  r})) 

(21) 

Rn  5  ((x,y)!(rn.j  ♦  rfl)  s  2R  •  2L} 

Therefore,  the  nodal  domain  cover  ft*  of  finite  element  ft*  Is  defined 
by 

fte»ftJU^+1  (22) 

where  (Fig.  A.l) 

a*  h  ((x.y)|r#  <  R  ^  (r^  ♦  Vl)/2} 

(23) 

'=  {lx.y)|(r#  ♦  r#+1)/2  <  R  <  r^) 
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Integration  of  the  governing  flow  equation  on  each  a*  Involves  the 
definition  and  Integration  of  the  thermal  conductivity,  K.  An  approach  to 
handling  the  nonlinearity  problem  Is  to  approximately  linearize  the  governing 
flow  equation  by  assuming  the  various  parameters  to  be  constant  for  small 
tlmesteps.  At.  For  K  set  to  a  quasi -constant  value  K*  In  ne  during  tlmestep  At, 
(19)  may  be  rewritten  for  the  nodal  domain  n*  contribution  to  subdomain  Rj 


» 

.  30 

' 

1 

30 

RKe  — 

ds 

,  9  < 

RS  —  dA 

. 

3n 

■ 

3t 

v® 

p®  n  r® 

a 

m  J 

.  JeS. 


(24) 


a* 


11  ‘j  “J 

Since  the  governing  heat  flow  equation  Is  radially  symmetric,  (24)  simplifies 
to  (see  Fig.  A. 2) 


o 

30 

RKe  — 

.  =  . 

RS  —  dR 

3R 

3t 

,e 

peOre 

e  1 

J  J 


.  JeSa 


(25) 


This  integrated  relationship  will  be  used  to  develop  a  subdomain  integration 
model  in  the  following  section. 


A. 3  NUMERICAL  SOLUTIONS 

In  the  derivation  of  the  finite  element  integration  statement  of  (25) 
for  Qe,  no  specification  of  the  character  of  the  state  variable  is  assumed. 

In  the  following,  the  state  variable  0  is  assumed  to  be  adequately  approximated 
by  linear  trial  function  0e  in  each  finite  element  fle.  Therefore  it  is  assumed 
that  0  =  0e  =  EL.  in  each  ne  where  the  L.  are  the  usual  linear  local 

J  J  J 

coordinates  in  0*;  and  0*  are  nodal  point  values  of  the  temperature  trial 
function  estimate  in  ne.  Oue  to  the  linear  definition  of  $e  in  fte,  all  scatial 


gradients  of  d*  are  constant.  Consequently,  several  well  known  domain 
numerical  solutions  of  (1)  In  n  embodied  In  the  finite  element  method  of 
weighted  residuals  result  In  similar  numerical  approximations  of  (1)  In  each 
Qe.  To  develop  these  domain  numerical  solutions,  the  following  description 
variable  Is  defined: 

34 

-  R$  —  (26) 

3t 

Although  the  Galerkln  method  of  weighted  residuals  used  to  solve  (1)  In 
each  fle  Is  well  known,  its  derivation  of  a  finite  element  matrix  system  is 
presented  in  order  to  develop  some  of  the  notation  and  simplifications  used 
In  the  subsequent  determinations  of  the  subdomain  Integration  and  Integrated 
finite  difference  analogs,  and  to  demonstrate  the  use  of  flux-type  boundary 
conditions  on  global  boundary  r. 


4  =  — 
3R 


34 
RK  — 
3R 


Galerkin  Method  of  Weighted  Residuals 
In  ne. 


4  Lj  dR  =  0 


n 


(27) 


generates  a  Galerkin  finite  element  matrix  system  for  approximation  of  (1) 
on  fie. 

Integrating  by  parts  reduces  (27)  to 


f  3b 

► 

34  dL.  34 

♦  Lj  dR  s  RK  —  L  , 

1  J  3R  J 

• 

RK  —  — ^  ♦  SR  —  L. 

3R  dR  3t  J 

dR 


(28) 


re  ne 


The  first  term  in  the  expansion  of  (28)  satisfies  flux  continuity  between 
finite  elements  and  Neimiann  boundary  conditions  on  global  boundary  r  in  a 
manner  similar  to  the  NOI  statement  of  (19). 


For  U.K)  =  (oe,  Ke)  in  ne  during  a  small  tlmestep  At  where  ♦*  is  an  assumed 
linear  trial  function  for  $  in  fle,  (28)  simplifies  to  the  Galerkln  finite 


element  statement 


A  dL1  3  A 

0  i  Ke -  R  — “•  dR  «•  S  —  R*eL,  dR 

3R  dR  at  ,  J 


where 


—  5  *  *e)/(Vl  ‘  r«> 


Integrating  (29)  determines  the  Galerkln  finite  element  matrix  system  for 
the  approximation  of  (1)  on  ne 


Se  Ke  ♦  Pe(2)ie  +  Qe(2)ie  =  {0} 


where 


r  +  r 
e  e+I 


1  -1 


-1  1 


e  St 

Pe  ( 2 )  =  -i- 
6 


Qe(2)  = 


s(te)2  ri  i 

12  1  3 


and  where  ie  =  (rg+j  -  r^);  and  ( tj>e ,  $e)  are  vectors  of  the  nodal  values 
time  derivative  o*  nodal  values  of  finite  element  e. 


Subdomain  Integration 

A  cover  of  finite  element  fle  Is  given  by  the  union  of  nodal  domains 
JeSe*  Th#  subdomain  Integration  version  of  the  weighted  residuals  process 
approximates  (1)  In  each  subdomain  Rj  by 

|  •  WjdR  =  0  (35) 


where 


1  .  R  e  Ri 


0  ,  otherwise 


♦  WmdR  ■  f  ♦  dR  ♦  f 


Thus  a  finite  element  matrix  system  is  generated  by  the  subdomain  integration 
method  for  finite  element  fie  by 


=  (0),  j  cS, 


From  (25), 


Using  (M)  =  Ue»  ♦*) 


* 

«  dR 

*  s  i 

► 

K* - 

R 

3 

-  s  — 

* 

R 4e  dR 

3R 

3t 

l  e 

p®_  r®  r® 

i 

.  je$.  (40) 


n 


j 


j  *r  “j  “j 


Integrating  (40)  gives  the  finite  element  statement  of  a  subdomain  Integration 
approximation  of  (1)  on  finite  element  0* 


f  $  +  !  <3)  !  +  9  (3>*  =  {0}  (41) 

where  Se  Is  given  by  the  Galerkln  element  matrix  of  (32) ,  and 

,e 


e  V 
Pe(3)  =  e 


Qe(3)  = 


8 

s(te)2 
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3  1 
1  3 

2 
2 


(42) 


1 
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(43) 


and  where  vectors  (<j>  ,<p  )  are  as  defined  previously. 

From  (40),  the  global  finite  element  matrix  system  determined  by  the 
appropriate  summation  of  each  ne  matrix  system  satisfies  Dirichlet  and  Neumann 
boundary  conditions  in  a  manner  similar  to  the  global  Galerkin  matrix  system. 

Additionally,  the  capability  of  representing  an  inhomogeneous  medium  by 
specifying  different  parameters  in  each  finite  element  fie  is  similar  to  the 
usual  Galerkin  finite  element  approach,  although  from  (40)  the  conduction 
parameters  are  necessarily  evaluated  at  the  midpoint  of  the  finite  element 
and  capacitance  parameters  need  to  represent  a  mean  value  in  the  control 
volume  associated  to  the  particular  nodal  point.  These  advantages  normally 
associated  to  the  Galerkln  finite  element  approach  can  also  be  developed  for 
an  Integrated  finite- difference  numerical  analog  for  the  approximation  of  (1) 


Integrated  Finite  Difference 


The  Integrated  finite-difference  approach  (Spalding,  1972)  can  be  extended 
to  the  solution  of  (1)  on  appropriately  defined  control  volumes.  The  usual 


control  volume  definition,  however.  Is  Identical  to  the  subdomain  definition 


cover  Rj  of  global  domain  n  given  by  (21).  Thus,  an  Integrated  finite- 
difference  approximation  for  solution  of  (1)  on  ne  where  the  trial  function  $e 
Is  assumed  to  be  linear  In  each  fle  Is  given  by 


rj-rjnr* 


■  (  S  —  R  (j>  dR  >  ,  j  e  S 


But  the  integrated  finite-difference  approach  equates 


♦j  dR 


Thus,  the  finite  element  matrix  system  Is  given  by 


se  *e  ♦  Pe(»)  *e  +  0e(«)  *e  =  to) 


where  S  is  once  more  given  by  (32);  and 


e  V  r 1  0 

Pe(»)  = 


2  0  1 


S(ie)2  ri  o 


0  (®)  i 


0  3 


Fro*  (40),  the  global  finite  element  matrix  system  determined  by  the 
appropriate  sumntlon  of  each  fle  matrix  system  of  (46)  satisfies  Neumann  and 
Olrlchlet  boundary  conditions  In  a  manner  generally  associated  with  the  Galerkln 
approach.  Additionally,  anisotropic  Inhomogeneous  mediums  are  similarly 
accommodated  as  with  a  Galerkln  or  a  subdomain  integration  analog  previously 
derl  ved. 

Due  to  the  similarity  of  the  three  numerical  approximations,  a  single 
element  matrix  system  for  the  approximation  of  (1)  on  fte  can  be  written  by 

Se  $e  ♦  Pe(n)  4*  ♦  Qe(n)  ie  -  (0}  (49) 

where  Se  is  given  by  (32);  and 


p  S(£e)2  r (n2  -  2n  +  3)  31 

0e(n)  = -  (51) 

4(2n2-n+3)  [_ <3n  -  3)  (3n2-3n+3)_ 

The  Galerkin  finite  element,  subdomain  integration,  and  integrated  finite- 
difference  numerical  analogs  of  (31),  (41),  and  (46)  are  given  by  (49)  for 
n  =  ( 2 , 3 ,*)  respectively. 

Extension  of  the  NOI  technique  to  cylindrical  coordinates  follows 
directly  from  (32),  (50)  and  (51).  The  extension  to  spherical  coordinates 

follows  from  the  tetrahedron  finite  element  determination  (Hromadka  and  Guymon, 
1983b)  and  the  above  results. 


A. 4  NOI  MOOEL  ANALYSIS 


The  previous  section  unified  several  numerical  techniques  Into  a 
single  finite  element  matrix  system  as  a  function  of  the  degree  of  mass 
lumping,  n.  The  question  remains  whether  an  optimum  n  factor  exists  such  that 
the  modeling  Integrated  relative  error  Is  a  minimum. 

In  this  work,  then  factor  developed  for  one- dimensional  diffusion 
problems  (Hromadka  and  Guymon,  1983b)  Is  used  to  test  the  NDI  model  accuracy 
for  radial  problems  where  analytical  solutions  exist.  This  technique  Is  based 
on  the  Fourier  series  expansion  of  a  particular  solution  to  the  one-dimensional 
diffusion  equation  In  a  small  homoaeneous  control  volume.  That  Is,  the  radial 
geometric  contributions  are  neglected  In  the  development  of  n. 

For  a  control  volume  Rj,  the  usual  process  of  normalization  reduces  the 
one-dimensional  diffusion  equation  to 

3J0  36 

-  *  —  ,  x  c  (0,1]  (52) 

3x2  3t 

where  0  is  a  normalized  variable  for  $,  and  variables  x,t  are  now  defined  as 
normalized  space  and  time.  It  Is  assumed  in  (52)  that  e(x*0)  *  e.  ,, 

J  * 

e(x-.5)  =  9 j ,  and  e( x  *  1 )  *  0j+j  where  are  normalized  nodal  values. 

Using  the  Crank-NIcolson  mldtimestep  advancement  procedure  to  approxi¬ 
mate  the  time  derivative,  the  nodal  equation  for  solution  of  0j  is 

S  -  “i*1  *  •j:5]  *  W-i  -  *•]  *  ej*i)  J  - 

(53) 

St  C(ej-J  •  9J->)  *  2nJ  (*r  -  9j)  *  -  9m)] 

J 


where  the  normalized  length,  |  i R j  1 1  -  7  i  1  Is  the  tlmestep  nwber;  and  nj  Is 

constant  during  normalized  tlmestep  At.  Equation  (53)  evaluates  all  modeled 

flux  tens  at  the  mldtlmestep.  For  other  time  derivative  approximations  such 

as  forward  or  backward  step  differencing,  a  similar  difference  statement  can 

be  developed.  A  solution  of  (52)  using  (53)  at  the  nridtlmestep  Is 

-  If  1  .-i- 

9(x,c)  «  -  2  9j_1  -  20 J  +  0j+J  sin  irxe 


+  l9j+l  *  9j-lJ  x  +  9j-l 

where  e  Is  normalized  time  measured  from  the  mldtlmestep;  and  where 
e  *  j  ©j  If  It  Is  assumed  that  all  effects  of  a  moving  boundary  value 

at  the  endpoints  Is  equivalent  to  holding  0  constant  at  the  mldtlmestep  boundary 
values,  then  (54)  represents  an  exact  solution  to  the  assumed  boundary  value 
problem. 

Holding  the  boundary  values  of  0  constant  at  the  mldtlmestep  allows  a 
simplification  of  the  NOI  nodal  equation  to 
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Solving  for  9j  and  9j+1  gives 
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Combining  (56a,b)  and  (55)  gives  nj  as  a  function  of  the  model  t lees tap  size 
by 

4ot  (l  +  e‘wlAt) 

Oj(At)  ■  ■  lii  ’  (57) 

j  -  4At(l*e  *  At) 

where  the  normalized  tlmestep  At  Is  ralated  to  the  global  model  tlmestep  At* 
by 


At* 


At 


^IIRjII2 


(58) 


where  0  Is  tha  mean  dlffuslvlty  ( D  ■  K/S)  for  Rj. 

From  (57),  the  mass  lumping  factor  lies  within  the  range 


8/ (w 2  -8)  <  Hj(At)  <  •  (59) 

and  Is  seen  to  be  a  function  of  tlmestep  and  element  size. 

To  test  the  success  of  the  n(At)  selection  technique,  several  radially 
symmetric  heat  transfer  (diffusion)  problems  where  analytic  solutions  are  known 
were  modeled.  Additionally,  the  derived  n  factors  of  2,3,  and  •  were  also 
tested  for  comparison  purposes.  The  measure  of  accuracy  used  Is  a  form  of  the 
l2  norm  of  the  error  given  by 


where  $  Is  the  test  problem  solution;  $  Is  the  approximation  value;  and  E  is 
the  error  of  approximation. 


In  this  stu4y,  six  different  boundary  value  problems  of  the  heat  flow 
equation  were  tested  with  various  values  of  the  tleestep  { Crank- N1 col son 
nethod)  and  finite  element  size.  For  each  test,  nodal  values  are  reset  to  the 
exact  nodal  values  after  each  tines tep  advancenent  In  order  to  better  test  the 
approximation  error  In  satisfying  the  flow  equation  rather  than  measuring  the 
accumulation  of  approximation  error.  After  each  times tep,  the  E  error  Is 
evaluated  and  stored  for  the  four  considered  n  factor  approaches,  and  the  factor 
which  resulted  In  the  minimus  E  value  (a  success)  Is  noted.  Consequently,  more 
than  150  test  problems  (5  times taps  and  5  element  spaclngs  for  each  boundary 
value  problem)  resulted  In  an  excess  of  20,000  times tep  advancements.  By 
dividing  the  number  of  successes  by  the  total  number  of  tlmestep  advancements, 
a  probability  of  success  for  a  n  factor  Is  estimated.  Figure  3  shows  the  success 
probability  for  the  n(At)  approach  as  a  function  of  normalized  tlmestep  size  and 
element  size. 

A. 5  CONCLUSIONS 

1.  A  unifying  numerical  model  can  be  developed  for  radially  symnetrlc  heat 
conduction  problems.  The  unifying  model  Is  based  on  the  straight  forward 
nodal  domain  integration  method.  The  resulting  model  Is  found  to  have 
the  capability  of  representing  the  Galerkin  finite  element,  subdomain 
integration,  and  Integrated  finite  difference  methods  by  the  specification 
of  a  single  mass  matrix  lumping  factor,  n. 

2.  The  global  matrix  system  composed  of  the  sum  of  all  NDI  elements  accom¬ 
modates  Dlrlchlet,  Neumann,  and  mixed  boundary  conditions  without  the  need 
for  special  finite  differencing  equations. 

3.  An  infinity  of  possible  domain  numerical  methods  are  possible,  and  can  be 
represented  by  the  NOI  model  for  specific  values  of  n. 


A  conputer  coda  based  on  the  Galerkln  finite  eleaent  Method  can  easily  be 
Modified  to  allow  a  variable  Mass  luMped  Matrix  systeo  and,  consequently, 
represent  an  Integrated  finite  difference,  subdoeeln  Integration,  and  an 
Infinity  of  other  doanln  Methods. 

An  Inp roved  ness  1 unping  factor  exists  (as  a  function  of  tlnestep  and 
finite  eleaent  size)  which  nlnlnlzes  approxlnetlon  error  nore  often  than 
any  of  the  other  considered  donaln  Methods.  The  probability  of  the  pro¬ 
posed  opt 1 nun  ness  leaping  systan  being  the  best  nunerlcal  Method  Is 
approxlnately  70  percent  for  the  no  null  zed  tlnestep  sizes  considered. 

The  Inp roved  Method  Is  developed  based  on  a  linear  trial  function  nodal 
and  a  Crank  Nlcolson  tine  advancenent  approxlnetlon.  Although  only  the 
radially  syametrlc  problen  Is  developed  the  extension  of  the  approach  to 
cylindrical  and  spherical  coordinate  problem  Is  straightforward. 
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